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ABSTRACT 



The basic outline of a digital computer (CRC-102A) program for the 
numerical solution of the approximation problem in the time domain is 
presented. The method used, in effect, is a real time convolution and 
Dirichlet series representation of the Laplace transform to effect a 
time-to-frequency transformation. The major portions of the digital 
computer program and the complete flow chart to solve the problem are 
included. 

The composition of a basic library of programs for the solution of 
all network synthesis problems and the integration of this program into 
such a library is discussed. 

The program for the solution of the roots of an n— degree poly- 
nominal and the Dirichlet series expansion program were completed during 
the author’s industrial tour at the National Cash Register Company, 
Electronics Division, Hawthorne, California. I take this opportunity to 
thank Mr. Henry Kent and Mr. Philip Sunday of that organization for their 
help in the writing and the debugging of the programs completed during 
my industrial tour* 

Special thanks are also due to Professor E. J. Stewart for his 
recommendations and constructive criticism with regard to the programming 
of the overall problem, and to Professors M. L. Cotton and J. B. Turner, 
Jr., for their guidance and help throughout the preparation of this 
paper. 



ii 



TABLE OF CONTENTS 



Section Title Page 

1. Introduction 1 

2. Statement of the Problem 3 

3. The Approximation Problem in the Time Domain 4 

4. The Network Synthesis Library 23 

3. Conclusions 28 

6. Bibliography 29 



iii 



LIST OF ILLUSTRATIONS 



Figure 

1. A Trivial. Case 

2. Derivation of h(t) 

3. A Triangular Inpulse Response 

4. Demonstration Problem 

5a. Three - Pole Approximation to h(t) -- Program 
Solut ion 

5b. Four - Pole Approximation to h(t) -- Program 
Solut ion 

6. The Basic Network Synthesis Library 



Page 

5 

8 

11 

19 

22 

23 

27 



iv 



1. Introduction 



The problem of network synthesis is a difficult and complex one. 
Consider the synthesis process to be divided into the two separate sub- 
problems of approximation and realization. There are many techniques 
for the attainment of the first goal, Z(p), the transfer function.^ 
Similarly, having attained Z(p), there are many procedures available to 
realize Z(p) in any one of several different circuit configurations. 

From the point of view of the electronics research engineer, a tool which 
would rapidly and accurately compute the circuit components of a network 
to provide a specified response, at a specified impedance level, in a 
specified circuit configuration would free the researcher from the tedium 
of designing the network and also from the temptation of accepting a 
•'second-best 11 solution. The designer would feel no compulsion to term- 
inate a problem rather than face another "messy" synthesis problem. 

The modern high-speed digital computer offers the engineer such a design 
tool, if an appropriate library of synthesis programs is available. 

Possibly the greatest advantage the computer avails its users is 
that once a certain type of problem has been solved, the method of solu- 
tion may be permanently recorded on punched cards or on magnetic tape or 
paper tape. The engineer now need only recall that such a problem has 
been programmed and refer to the catalogue of computer programs rather 
than to a textbook to refresh his mind on how to solve the problem by 
paper and pencil. 

*For a comprehensive summary of approximation methods, the reader 
is referred to Trans, IRE, Prof Grp on Circuit Theory, Vol. CT-1, No. 3, 
of Sep 1954. In addition to the summary, the first paper presented by 
Stanley Winkler, has a definitive bibliography of 240 items on network 
synthesis . 
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It is the purpose of this paper to present a program for the solu- 
tion of the approximation problem in the time domain and to show how this 
program may be integrated into a library of programs which use convention 
al methods for the solution of the more common network synthesis problems 
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2. Statement of the Problem 



The approximation problem as considered in this paper is felt to be 
defined as: 

Given the excitation available or expected at the input terminals 
of a passive four-terminal network, and the response or output required, 
or the impulse response, where these parameters are expressed in graphical, 
or general analytic form; find a rational function expressed as two poly- 
nominals, where one is understood to be the numerator and the other the 
denominator, which describes the transfer function within the limitations 
on the allowable deviation of the output or on the number of elements 
to be used (but not both), specified. Provisions are to be included for 
compensation of dissipative effects. The polynominals will be expressed 
either as polynominals or as the quadratic factors of polynominals. The 
function is to be physically realizable. 
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3. The Approximation Problem In the Time Domain 

Heretofore, approximation problems in the time domain have Involved 

lengthy and inaccurate frequency domain expansions which do not always 

lead to realizable system functions. These methods Involve complicated 

1(3) 

integrations, or at best, Taylor series expansions. At any rate, 

they are not suitable for the flexible program requirements we desire. 

What we would like to have is a method which would rapidly transform the 

input 'output characteristics of one network into a realizable system 

function within any desired degree of accuracy, or alternatively, an 

optimum realization for a specified number of elements. The numerical 
1(3) 

method of Ba Hll is tailor-made for these requirements. Furthermore, 
not only will it operate on input-output relationships, but also there 
is no reason why an Impulse response which has been obtained analyt- 
ically cannot be plotted and introduced into the program at the appropri- 
ate point to produce the transfer function for all types of networks, 
i.e., predictor networks, smoothing filters, interpolation filters, com- 
pensation networks and so on. In fact, we may say if a physically realiz- 
able h(t) can be defined, we can find the system function. The designer 
must have a plot or sequence of ordinates which describe the h(t), or 
the f^(t) and f Q (t), and use this plot or sequence to enter his problem 
into the approximation program. 

Briefly stated, Ba Hli's method is to find, by a process of synthetic 
division, a sequence ( £hj ) which represents the Impulse response of 
the network which in turn relates the input waveform to the desired output; 

^Freddy Ba Hli, Network Synthesis by Impulse Response for Specified 
Input and Output in the Time Domain, Tech. Rpt. No. 261, MIT Res. Lab. 
of Electronics, July 31, 1953 (hereinafter abbreviated as NSIR-BH). 



4 



and then to calculate the psuedo-system function H*(s), from £h] by a 



simple power series expansion. The H*(s) is a polynorainal in s which is 
P m t s ) 

the ratio of the which is normally considered to be the expression 

of the system function. In effect, we convert the convolution integral 
into its component operations and use an iterative substitution method. 
We shall not reproduce the theoretical reasoning and proof which was so 
neatly done by Ba Hli, but we shall demonstrate the ease and simplicity 
of the method with a numerical example. 

Before proceeding with a complete example, let us show how the 
synthetic division process yields the impulse response in a couple of 
trivial cases, the results of which are well-known. Consider Figure 1, 




Step 

i. i. i. i. 



Ramp 

. 2 , . 4 , . 6 , 



1 . 



• 8 , 1 , 0 , 



A Trivial Case 
Figure 1 
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the input f,(t) we take as the unit step whose Laplace transform is — ; 

L 8 

f (t), the output, is the unit ramp. Its transform, of course, is ^2. 
o s 

If we sample the input and output at intervals of .2 of a time unit, the 

sampled ordinates for f^(t) are ^1 , 1, 1 , 1 , 1, 1 ^ , and for f Q (t) 

|.2 , .4 , .6, .8, 1.0, 1.4 , 1. 6 Now dividing ^f Q (t)^/ ^ in the manner 



illustrated: 1,1,1, 1,1,1..../. 2,. 4, .6, .8, 1.0, 1.2, 1.4 

we have £hj^ « .2, .2, .2, .2, .2 and if we divide ^hj by At, we pro- 
duce the sequence 1 , 1, 1, 1 , 1 ,... .again the unit step. This, of course, is 
completely analogous to the process in Laplace transform operations, viz: 



or H(s) - W/ 





m 1/s 



We can apply the process in reverse, for example, |f^(t)l x {V} ^ 
■ ^f Q (t)^, we use un * t doublet as ^ we should differentiate 
the output. Consider the sequence below which defines a sine wave, with 
At® .1 

{f^ - .1 .199 .296 .389 .479 .564 .644 .717 



1 1.99 2.96 3.89 4.79 5.64 6.44 7.17 

-1 -1.99 -2.96 -3.89 -4.79 5.64 6.44 

[fji 1 .99 .97 .93 .90 .83 .80 .73 

jf J (Cont'd) .783 .841 .891 .932 .964 .985 .997 



7.83 8.41 8.91 9.32 9.64 9.85 9.97 

7.17 7.83 8.41 8.91 9.32 9.64 9.85 

(Cont'd) .66 .58 .50 .41 .32 .21 .12 
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and is a sequence which defines the cosine which is as expected 

since the unit doublet acts as a differentiator. Thus we demonstrate 

the ability of this ^h^ to perform as an operator equivalent to H(s) 

in the frequency domain on functions defined as sequences in the time 

domain. The transformation of £h^ to H*(s) follows logically from 

the definition of the convolution integral: 
t 

f 0 ( 0 - / f t (t) h (t -T) d T 
( 3 ) 

Ba Hli has shown that H*(s) can be represented by the Dirichlet 
power series, 

oO 

H*(s) ■ 5 a exp (-T s) 

L. n n 

n«l 



where the are the areas in our £h| 



^ sequence, and the are de- 



fined by 



T 

n 



n A t - 



A t 



that is, the T n are measured from the origin to the center of the sampl 
ing interval under consideration at the moment. 

We may approximate H*(s) as closely as we please by taking our 
sampling points at smaller and smaller intervals. 

To demonstrate the complete process and provide a basis for discus 
ing some of the more obscure points, let us take another example and 
follow it through more closely. Figure 2 shows the f^(t), the desired 
f Q (t) and the resulting h(t) required, (& t « 0.1): 

{f^tij-l, 1 , 1 , 1 , 1 , 0, 0 , 0 , 

{f o (t)J-.l, .2, .3, .4, .5, .4, .3, .1, 0, 0, 0, 

^ J *•!* •!» • 1 * • 1 » • 1 > 0 , 0 , 0 , 

In tabular form, we have 
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I 



Tn 

.05 

.15 

.25 

.35 

.45 



a 

n 

.1 

.1 

.1 

.1 

.1 



Substituting into our expression for H*(s), we have 



H*(s) » .1 [l - .05s + 

i H ic L • l! 
+ .1 j 1 * . 15s + 2 • 



( .05s )' 
31 

(.15s)' 

3’ 



+ .1 [1 - .25s + . Ujjl 



3! 



+ a . .35s + 



+ .1 [1 - .45s + 



(.45s)' 

3! 



(♦05s) 

4! 

Ldp j j! 

4: 



L-25-S , ). 

4 : 



(*35s) 
4 1 

( ,45s 
4.' 



- -D 

- ...] 

- -3 



« .5 - ,125s + .020625s - .00255208s + .0002517968s - 

If we compare this result with that given by Laplace transform methods, 
we have 

^[ f i<t)| - J [l * exp (-l/2s)] 

^l) 0 (ti] - j 2 Cl ‘ exp (-l/2s)] 2 

I[h(t)J “^^° <t ^[f 1 (t)]“ i[l - exp (-l/2s] 

and the power series expansion for ^[h(t)j is 

.5 - .125s + .020833s 2 - .002604166s 3 + .0002605166s 4 
which shows the first two terms by Ba Hli's method agree exactly, a 1Z 
error in the third term, 2% in the fourth, and about 4% in the fifth. 
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Before proceeding let us see how this accuracy may be improved. 
Clearly increasing the number of terms taken and maintaining the same 
interval will have no effect on the formation of the second, third, and 
fourth degree terms whose accuracy we wish to improve and intutively one 
feels that increasing the number of sample points will; therefore let us 
increase the number of sample points to, say, 8. 

With f (t), f (t) and h(t) as before we tabulate, for At » .0625: 



Tn 


a 

n 


.03125 


.0625 


.09375 


.0625 


.15625 


.0625 


.21875 


.0625 


.28125 


.0625 


.34375 


.0625 


.40625 


.0625 


.46875 


.0625 


« .0625 


[l - .03125. ♦ . ( > 93125 . ) 3 + (,.03 125.)* . 


.0625 


[l - .09375. ♦ <- 09 » 5s > 2 - (-°”; 5s > 3 - LmnsA . 



■i 



and so on; the resulting expression is: 



H*(s) - .5 - ,125a + .020635719s 2 - .002631105s 3 + .000257032861s 4 - ... 
and the new expansion for 8 sample points has produced errors of 0% in 
the first two terms, .95% in the third term, 1% in the fourth and 1.3% 
in the fifth. Therefore, the increase in the number of points has 
obviously resulted in an overall improvement. A further increase to 
10 points, produces accuracies of 0%, 07., .25%, .5%, .8% and 1.25%. 

The foregoing results are tabulated in Table 1 for the purpose of 
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easy comparison. 



We recognize that the errors in the higher order terms are caused 
by replacing the Laplace-Stiel jes integral^ by a truncated series summa- 
tion. This is equivalent to saying that we have not computed the true 
"moments" of a^ about the origin. 

Consider Figure 3, we have here a triangular impulse response. Let 

us calculate the second and third moments of this curve and compare the 

results with the power series expansion of the Laplace transform and 

also the Dirichlet series expansion. We let b^ be the second ‘‘moment", 

(it would be more precise to call this by a different terra since it is 
M 

actually 1 . where M is the moment), b equals the third; thus, 

j! j 3 



h(t) 




0 ^ t ^ 1/2 
1 / 2 ^ 1 
elsewhere 



A Triangular Impulse Response 
Figure 3 



l NSIR-BH (3) , pg. 15. 
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True Expansion .5 -.125s + ,0208333s - ,002604166s + ,0002604166s - ,000021701388s 
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1_ 1 


r 4| 1/2 
fx_ 


3 

X 


1 


4 

X 


- IT 1 


L4 '0 


+ 3 


1/2 


‘ 4 


“ 2 ( 


1 1 
64 + 3 


1 

" 24 


1 

' 4 + 


k> 



.03645833 



b 3 = jy {_ x 4 dx + n 2 ( (1-x) x 3 dx] 



I_ r 2L 

3 ! L5 



I 1 ' 2 < 

x 

'0 + 4 



1 

1/2 



- “ 1 

5 1/2J 



I ,_L_ l i_ l l v 
“ 6 *•160 + 4 “64 " 5 + 160' 

- .00781250 

The trancendental Laplace transform for this function is: 

H (s) - “7 [l * exp (-l/2s 2 2 

whose power series expansion is 
H(s ) • .250 - 0.1250s + .03645833s 2 - .00781250s 3 
The Dirichlet series results from 



'n 

.1 

.3 

.5 

.7 

.9 

H(s) » £ a n 
n-1 



-T.« 



n 

.02 

.06 

.09 

.06 

.02 



Therefore, considering only the b£ and b^ terms, we compute: 



13 



I 



2 2 
1 3 

b 2 « .02 x + *°6 x + 


. 5 * 

.09 x ~T + .06 x 




- j (.0002 + .0054 + .0225 


+ .0294 + .0162) 




« .03685 






Similarly, is computed as 

b 3 » .008008 






Therefore: 






"Moment" 


Laplace 


Ba Hli Approx. 


b 2 .03645833 


.03645830 


.036850... 


b 3 .00781250 


.00781250 


.008008. . . 



The obvious method of reducing the error is to take a smaller in- 
terval and thereby approach the true integration to a closer degree and 
since this will be a one time only calculation, we are not adamant about 
making extra computations to improve the accuracy • We also notice that 
by making our intervals smaller we reduce those errors which may arise 
from the process in which we obtain the £hjj For a more complete dis- 
cussion of the reduction of the inherent error in Ba Hli*s method see 
Appendix II. 

Having obtained a psuedo-sys tem function, it yet remains to deter- 
mine its equivalent as the quotient of two polynominals , the degree of 
each as yet unknown. We equate our psuedo-function, H*(s) to an express- 
ion of the form 



H*(s) 



P(s) 

Q(s) 



p 0 + p l s + P m s 



m 



q 0 + v + 



n-1 n 

q ,s + q s 
n n-l n n 



and solve for the unknown p*s and q's, since these coefficients express 



a ratio we may arbitrarily set any one of them we choose, (normally q^) 
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equal to one without any loss of generality. 

We notice that the series of equations which results from this opera- 
tion will always have the form: 

p o ■ Vo 

p i - Vo * Vi 

p 2 ■ Vo * Vl + V2 

0 . h 3 q Q + h 2<ll + * h 0 , 3 

0 . b 4 q 0 + hjqj + h 2 q 2 + * h # q 4 

0 - h 5 q 0 + h 4 , 2 + h 3 , 2 * h 2 , 3 ♦ h t q 4 

0 - ty*o * h 5 q 2 * h 4 q 2 + h^., 4 h 2 q 4 

For m « 2, n * 4; r, the relative degree between P(s) and Q(s), is 
obviously two. We see that we have seven equations in eight unknowns, 
which is as it should be for the expression of a ratio, and hence, re- 
quiring that we set q^ (or any other coefficient) equal to some conveni- 
ent value. The solution of this portion of the problem is considered in 
detail in Appendix III. 

We choose r, the relative degree of numerator and denominator by 
observing the behavior of h(t) around 0. Since we know from the initial 
value theorem of Laplace Transform Theory,** 
lim h(t) « sH(s ) 
t — >0 s — 

we may say that the function must behave as 
m 

as t approaches zero. 



P m s 

m 



V 



m 



Si. F. Gardner and J. L. Barnes page 265. 
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Thus , if h(t ) starts out as a step 

p/ Q \ / 

v . — * 1/s as s approaches infinity; or is h(t) starts 

/ W v s ) 

out as a ramp as S a PP roaclies infinity. 



To program the computer to sense this limit and thus decide whether h(t) 
is an impulse, step, ramp, or higher order function, we first determine 
h^, i.e., h(0). Since in some cases this may result in division by zero, 
we perform a first-order backward interpolation, using h^ and h^ to find 
tig. This is approximate but certainly is sufficiently accurate to deter- 
mine whether or not h(t) starts as a Jump function or as a ramp, which is 
the primary decision in the "selection of relative degree routine" (see 
Appendix I, page I-l). If it is decided at this point that h(t) is a 

first-order approximation to a ramp function, we must further decide 

3 4 

whether it is truly a ramp or a higher order function, i.e., l/s , 1/s , 
etc. These decisions are based on the fact that the n— derivative of an 
n-order polynominal is a constant, the process is flow charted for 
+1 ^ r ^ -4 on page 1-2-1 of Appendix I. The routine may be easily ex- 
tended to encompass all possibilities which may arise. 

We know that we can obtain any desired degree of accuracy in the 
terms of our power series expansion and further that we can correctly 
determine the relative degree between our numerator and denominator poly- 
nominals. We have shown that we can compute the coefficients of our H*(s) 
and we know that since this is a ratio with arbitrary coefficients we may 
multiply the expression by the impedance level factor to set any required 
impedance level. But what of the approximation error, i.e., the error 
caused by not taking an infinite number of terms to approximate tran- 
cendental functions ? 



Ba Hli has shown that in the case where the Laplace transform of 
h(t) is rational, it will be recovered exactly. NSIR-BH,(3) pgs. 41-44. 
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N o 5 *1 + 1.5(5) + 2 
Therefore, we can write 



13.5, which we take as 14. 



Pq + PjS ♦ P 2 s 2 + ... + P 12 s 12 + p 13 s 13 11 27 

2 13 14 “ h 0 + h l S + *•’ h ll S + h 27 S 

q 0 + q i S + q 2 S + + q i3 s + q i4 s 

this expression leads to 23 equations in 24 unknowns, one of which, custo- 
marily q^, we set equal to 1, without error. To avoid filling sheets of 
paper with useless symbology, we shall terminate this example here and 
turn to another pair of functions. 

We shall compute, using the program, approximate impedance functions 
for a pair of input-output time functions and check the computer solution 
analytically. To do this, we choose the functions illustrated in Figure 
4. 

f^(t) is a rectangular pulse 
f i (t) * 1 0 £ t — 1 

0 elsewhere 



f^(t) is an inverted sectioned parabola 





f 0 ( £ ) - . 


5t 2 




0 - t - 1 










m - 


t 2 + 2.0t - 


1.5 


1 & t ± 2 










■ 


5 (t-3) 2 




2 - t — 3 










O 0 






elsewhere 








for At 


* .1, the following sequences obtain 








£i<‘>3 


- i, i, 


1. 1. 


1. 


1. 1. 


1, 


1. 


1, 




- .005, .020, 


.045, .080, 


.125, 


.180, .245, 


.320, 


.405, 


.500, 




.590, .660, 


.710, .740, 


.750, 


.740, .710, 


.660, 


.590, 


.500, 




.405, .320, 


.245, .180, 


.125, 


.080, .045, 


.020, 


.005, 


.000, 



Synthetic division of -jf Q (t ) / f ^ ( t ^ yields : 



18 




tti:L functions 



Senicaatrat iBtL 



i 


tiTj; 


■ ' r ■ - 

ir" _ : 


i i-L. 

'M: 
-1-1-1 . 


“ X’ 
444- ■■ 

■;4 


i i i . 


-t tj.. . 


W : - 


w. 


:±i . 


••■4 ♦ - 

_0 y _ 

..jJXl . _ 

fH 




J_! 1 _ 

U! ' 


’ ri-f " 



19 



- .005, .015, 


.025, 


.035, 


.045, .055, .065, 


.075, .085, .095, 


.095, .085, 


.075, 


.065, 


.055, .045, .035, 


.025, .015, .005 



w 



which describes h(t) as the triangular function also shown in Figure 4. 
We know that the Laplace transform from f^(t) is: 






i- /i " s \ 

i < i - e > 



and the Laplace transform for f (t) is 



o 

s -.3 



£&.<*>] - 



Thus H(s), which should be £f o (t)J / jjf ^ ( t )] or — J , which is 

2 

the transform of a triangular pulse . 

We expect H*(s) as computed by our program to be a good approximation 

-s-, 2 



to the power series expansion of 

-S-. 2 



M • 



L^J 



1.000000000 - 1.000000000s + .583333333s 

- .250000000s 3 - .086111111s 4 - .025000000s 5 
+ .006299603s 6 - .001405423s 7 + . . . . 



- 1.000000975s - .583340832s 2 - .250007019s 3 + 

- .025001670s 5 - .006300172s 6 - .001405587s 7 + 



The program approximation with single interpolation, i.e., At/2, produced: 
H*(s ) o .999999999 
.086115077s^ 

If we now solve for the impedance function using n » 3, and r =» -2 
(since h(t) is a ramp function in the vicinity of zero); 

0.999999998 - 0.049971105s 



Z(P) 



1.0 + 0.950029868s + 0.366689961s 2 + 0.062506119s 3 

» ( - .049971105 ) (s - 20.011564640) 

(s + 2.277182620) [(s + 1.79464120) 2 + 1 .950589991 2 ] 

^Gardner and Barnes ^ 6 \ Appendix A, No. 6.08. 

Gardner and Barnes' 1 , Appendix A, No. 6.07. 



or; 
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The inverse Laplace transform of this function is: 



4 .4132006227 e' 2 - 27718262t + 4.465570186 e' 1 • 79464 120t sin (l.95058999t-81 .217°) 
The plot of this inverse superimposed on the desired h(t) is shown in Figure 
5a. 



If we desire to improve this approximation we go to a 4-pole network. 
The program computes the system function as: 

2.67414555s 2 - 8.1622705s + 135.9105227 

s 4 + 10.59828491s 3 + 51.1404969s 2 + 127.7483847s + 135.9105230 



which we expand as : 
1.0681279249 



fr 



fs + 23.98432161) 



(s + 20.00984647) 



l [(s+3. 017294581 ) 2 +(. 932721456 ) 2 J [(s+2 .281847873 ) 2 +(2 . 901655110) 2 ] 

and determine the inverse Laplace transform to be: 






24.010971222 e " 3 .01729458U gin ( >932 721456t + 2.547°) 

-6.612687409 e _2 * 2818A7873t sin (2 .901655110t+ 9.295°) 

This inverse is plotted in Figure 5b. Note the improvement with the 
addition of the other pole. The approximation can, of course, be further 
improved by adding still more poles. 

Suffice it here to say, that the existing program will solve for the 
impedance function, given the time-series data describing the input and 
output waveforms. 

Before concluding this section, we estimate the required solution 
time for the total approximation problem in the time domain. 

We assume the following conditions obtain as a basis for our estimate: 
1. h(t) is not a closed function, i.e., K » 4 (see Appendix II). If 
h(t) is closed, that is, if h(t) is completely described in 4 or less time 
units the running time for the synthetic division routine and for the cora- 
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putation of the psuedo-system function, H*(s), will be proportionately 
reduced. 



2. H(s) starts as a step function, i.e., r » -1. This is the worst 
case we encounter as regards computation time for the psuedo-system funct- 
ion, since we must compute 2<n + 1) + r terms of the expansion of H*(s) 

in order to solve the system of simultaneous linear equations. 

3. The interpolation decision is such that one intermediate inter- 
polation calculation is performed. This is the intermediate condition, 
and requires twice the computation time for a no-interpolation decision 
and one-half the time for the four-point interpolation. 

We consider the overall program (for time domain approximation) to 
be sub-divided into three, (four if Z(p) is desired in factored form or 



if compensation is to be provided) major computations. Thus we have.- 



Sub-division CRC time 

Synthetic Division to 5 m i n * 

produce ^hj 



Word-times 



72 x 10 



4* 



Dirichlet power series 
expansion H*(s) 

Matrix Inversion to 
compute Z'p 

Factor Z'p to apply 
loss compensation 



n x 3 min 
.064 n^ min 
(2n-l) x 9.4 rain 



(n-1) x 44 x 10* 



.96 n^ x IQ** 



(2n-l) x 144 x 10* 



♦Estimated, all other values measured. 

For a 10-pole network, this gives a solution time of 278 minutes; for 
a 15-pole, 538 minutes; for a 20-pole, 945 minutes. 
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4. The Network Synthesis Library 

Figure 6 shows the block diagram of a network synthesis library for 
the automatic computation of the approximation functions for a variety of 
filters. This diagram is by no means considered to be all-inclusive, it 
is only intended to show how such a library would be assembled. 

It should be obvious from the diagram how our various blocks are used 
together to provide a large number of filter characteristics from a rela- 
tively small number of programs. All blocks would have common input, 
output and carry* locations, and programs represented by the blocks would 
have bootstrapping links to call in the next block in order to proceed 
with the computation. 

In addition to the "f requency" filters shown, we would consider it 
desirable to provide in the basic library: 

a. Tschebycheff Stop-Band Ripple Filters 

b. Tschebycheff Equi-Ripple Filters 

With regard to the latter, Byron J. Bennett^ has developed an inter- 
esting iterative technique which seems to be much simpler than Alexander 
2 

J. Grosman f s method for avoiding the elliptic function. 

Note the provision for the realizibility check. It is considered that 

the inclusion of such a routine is mandatory in order to save computer 

3 

time. This routine would compute n, the number of elements required to 

*We use the term carry to characterize that data which is not used 
in the present computation but must be "carried" for use in later routines. 

^Byron J. Bennett, A Note on Circuit Synthesis, p. 61.^^ 

2 

A. J. Grosman, Synthesis of Tschebycheff Parameter Symmetrical 
Filters, Proc. IRE, Vol. 45, No. 4, pp. 454-474. (^) 

3 C 13 \ / 4 \ 

Using Grosman's' * equation (38) or see Darlington^ . 
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provide the required selectivity parameter (or discrimination parameters) 
and meet the ripple requirements , if n is specified and the requirements 
demand more elements than is allowed the computer should so indicate to 
the operator. It is estimated that this information would be available 
within not more than ten minutes from the start of computations, and 
therefore establish a doctrine that if such a program is used, the opera- 
tor must wait until this point in the program is passed before allowing 
the computer to proceed unattended, and by so doing be assured that the 
results will be realizable as a physical network, and that computation 
time will not be wasted. 
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BASIC NETWORK SYNTHESIS LIBRARY 
Figure 6 
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5. Conclusions 



The results of the tests on the routines programmed to date indicate 
that a library of network synthesis routines employing conventional tech- 
niques is practicable. 

The use of such a library should result in considerable savings in 
time and effort for those people concerned with network design. 

The accuracy obtainable is more than sufficient to guarantee excellent 
approximations (and undoubtedly synthesis) if the number of elements to 
be used is not severely restricted. 

The time required for solution is not excessive in comparison with 
the magnitude of the problem and the accuracy of the results. The employ- 
ment of the library of programs will require a working knowledge of the 
computer, unless magnetic tapes are made up from the “building blocks** 
to provide specific types of network designs. In this case, detailed 
instructions may suffice to permit persons not familiar with the computer 
to use the programs to good advantage. 

The installation of the CDC - MK 1 Computer will enhance the value 
of this library despite the conversion problem. The increased speed of 
computation and shorter programs which will be possible due to the “built- 
in** floating point arithmetic will make such a library even more valuable. 
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APPENDIX I 



It is the intention of this appendix to explain with a minimum of 
theoretical reasoning the overall flow chart of Ba Hli's method. 

The following inputs are required: 





i. f o 


<t)] - 


time sequence of output function 


(in octal 


or decimal). 


2 * J 


l f i 


<t)j - 


time sequence of input function 


(in octal 


or decimal). 


3. 


N d 


- 


number 


of 


discontinuities . 






4. 


N 

s 


- 


number 


of 


points of inflection. 






5. 


N 

e 


- 


number 


of 


elements. 






6. 


t 


- 


number 


of 


f^(t) ordinates used. 







7. Octal or Decimal input data indicator. 

Annex 1 shows a block diagram flow chart which considers only inputs 
and outputs to major program subdivisions, and major logical decisions 
required. 

We allow for either decimal or octal digit input sequences, in either 
case they are converted to binary floating point. If the number of ele- 
ments (N e ) is not specified, compute AS « 1.5 + N g (see page 17,raain 

text). 

The synthetic division process is self-explanatory, the quantity 
"t 1 *, the number of f^(t) ordinates, is required in this block. 

Annex 2 shows a detailed flow chart for the determination of ,, r H , 
the relative degree between numerator and denominator. The first order 
approximation used to compute h^ is considered to be sufficiently accurate 
to ascertain whether or not h(t) starts from a finite value, or from 
zero. The ,, 0 M used in the decision boxes is actually slightly greater 
than zero, about .025, to avoid any erroneous decision that h^ is greater 
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than zero, when it is only greater due to the error introduced by the 
first order approximation. Note that we have deviated from Ba Hli's 
method of determining the behavior of h(t) near t « 0 in that we compute 
derivatives . 

We can now determine, D^, the degree of the denominator, and N^, 
the degree of the numerator of Z(p) and the number of terms of H*(s), 

H^, we must compute to solve for Z(p). 

We can make a somewhat arbitrary decision as to how accurately we 
should compute our H(s), i.e., what subdivisions we should use in our 
"moment 11 computing scheme (see Appendix II), depending on the number of 
terms we shall have in our denominator. 

The next two appendices give the program details for the solution 
of H* ( s ) , and Z'(p). 

The solution of Z(p) is now complete, unless we desire to compensate 
for dissipative losses. This is easily accomplished after factoring the 
numerator and denominator polynorainals . We preshift our pole and zero 
pattern to the right an amount " <*•", unless this shift causes a pole or 
zero to move into the right half plane and thereby violate the conditions 
of physical realizability. We compute 
^ « R/2L 

2L/R - 1/c^ 

coL/K » *V2^ 

2Q /oj m cA. 

^ » Q/rrf 

and thus we compute < as a function of Q and f, where f is the lowest sig- 
nificant frequency involved and is entered by the designer as input data. 
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Realization Programs 



OVERALL FLOW CHART 
OF 



TIME DOMAIN NETWORK SYNTHESIS PROGRAM 



Selection of Relative Degree 
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APPENDIX II 



CALCULATION OF H*(s) FROM THE [h] SEQUENCE 



We are now required to form the power series expansion of: 



s -T. 

Z V 1 

i-I 

As was mentioned in the main body of this paper, replacing the Laplac 
Stieljes integral by the summation indicated above, resulted in an increas 
ing error in those terms higher than the first degree. We recognize that 
this error results from the failure to compute the exact moment of i*h{ ^ 
about the origin. 

The true moment of an area defined by f(x) is: 

°0 



M 



j 



j 



f(x) dx 



which we have replaced by the summation 

4 r 2 t ' h J 

Hj(app) . 2^ [-J 5 -J \ 



where A. * incremental area or h and S a NK (K * 1,2,3,... ) 

i a^ 

For our program, we limit K to 4; that is, we require the impulse 
response to be defined over the range 0 to 4 time units. The input and 
output functions may be sealed so this will be adequate for most problems. 

It might be well before proceeding to take a close look at the manner 
in which we form the coefficients for the Dirichlet power series expansion 
We may proceed in such a manner as to form one term at a time, that 
is, referring to Figure II-l, we may sum first down the j » 0 column, then 
the j * 1 column, etc. The operation count for this procedure is: 

N "a" for the first term, 
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N "a" + N ,, m n + "m" for the second (if we factor out the 1/2N until 

the sum is complete), 

N "a" + N(j + l) ,, m M + (j + l) M m n for the j*"* 1 term, 
a total of 

& 

(j + 1) N"a" + N+l 2 (i + 1) "m" 

i»l 

(j + 1) |n fa" + l^r + l) "»"] +( + l) "m"j 

where "a" stands for floating point addition and "m" stands for floating 
point multiplication. 

For N « 40, j « 20; this becomes 840 "a 11 + 11040 "," and if we assume 
that floating point addition requires twice the time that floating point 
multiplication does, a equivalent total of 12,720 floating point multi- 
plications . 

If, however, one proceeds out the rows and forms all the partial sums 
associated with A^, then , etc., for i « 1, j » 0, to j, the count per 
row is : 

"a", 2"m" + "a", "a" + "m", "a" + M m n , "a" + "m" and the 2N in 

the denominator is included. 

The count is the same for all rows, therefore, the total is: 

NU+IX'WW), which for N =« 40, j « 20, is 840 "a” + 840 "m", or 
based on our previous assumption, a total of 

2,520 floating point operations. 

There is, however, some additional "call-up" and "put-away" address compu- 
tation which degrades this efficiency somewhat. Nonetheless, it seems 
quite evident that a saving of at least 507. can be effected in this way. 

With this saving we can afford to reduce the error in the final co- 
efficients by introducing an interpolation routine which effectively 
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doubles or even quadruples the number of sample points used. 

The interpolation precedure makes a straight line approximation be- 
tween sample points, and computes the value the intermediate sample points 
would have. Therefore, it is able to use a smaller increment on the 
abscissa and thus approaches a true integration more closely than would 
be otherwise possible. The results of some runs using various M degrees n 
of interpolation are presented in Table II-l for comparison purposes. 

The flow chart and program are appended hereto. 
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N 



•i— A 

2N 1 



( 2N ) A i 



( 2N^ A 1 



( 2N > A i 



- — 

2N 2 



( 2N ) a 2 



— a 
2N 3 



( 2N ) a 3 



2N A 4 



— a 
2N 5 



^ 2N 



< 2N ' 



^ 2N ' 



21-1 ^ 

A i 



N 



( 



2Nliw 
2N )A n 



x ^ , 2N-i . . , , 2N-i v 2 . , «> ,2N-1. 3 , , , ^ , 2N-U j . 

H (s) - 2 A i + Z ( 2N * A i + 4 ^ 2N * A i + Y ^ 2N ^ A i ’ * * + ^ ^ 2N * A i 
i i i 1 l 



The factorials which appear in H(s) are yet to be supplied, and are 
brought in with the alternating sign according to whether j is odd or even. 



Formation of H*(s) 
Figure II-l 
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DIRICHLET POWER SERIES EXPANSION FOR PSEUDO-SYSTEM FUNCTION 



ADD 


INST 


M 1 


M 2 


M 3 


REMARKS 


0300 


05 


3000 


3000 


0550 


Load buffer 


0301 


3** 


3000 


2100 


2000 


and clear putaways 


0302 


17 


2010 


3000 


0305 


If TS 1 up, 0305 (l) 


0305 


30 


0556 


2100 


0545 


No, Set t to 1 


030 U 


34 


3000 


2100 


0311 


and skip to 0311 


0305 


17 


2020 


3000 


0310 


If TS 2 up, 0310 


0306 


30 


0541 


2100 


0545 


No, Set t to 2 


0307 


34 


3000 


2100 


0311 


and skip to 0311 


0310 


30 


0542 


2100 


0545 


Set t to 3 


0311 


35 


0536 


0545 


0535 


N exp + t to N' exp 


0312 


32 


0510 


0516 


0314 


Set h* pickup 


0313 


31 


0533 


0574 


0551 


"- 1 " to SA 


0314 


31 


0120 


0121 


0523 


h t to TF 


0315 


3** 


0545 


0532 


0321 


If t * 1, 0321 


0316 


31 


0523 


0564 


0530 


TF to TA 


0317 


30 


0533 


2100 


0547 


1 to H 


0320 


34 


3000 


2100 


0413 


Skip to 0413 


0321 


34 


0545 


0534 


0350 


If t * 2, 0550 


0322 


33 


0572 


0533 


0326 


No, if SA * "0”, 0326 


0323 


36 


0523 


0533 


0521 


No, TF/2 to A 


0324 


30 


0564 


2100 


0562 




0325 


34 


3000 


2100 


0333 


Skip to 0333 


0326 


31 


0523 


0564 


2000 


Load TF 


0327 


31 


0524 


0565 


2002 


Load TE 


0330 


34 


3000 


2100 


1516 


TF - TE 


0331 


36 


2000 


0534 


2000 


(TF - TE)/4 


0332 


31 


2000 


2001 


0521 


(TF - TE)/4 to A 


0333 


31 


0523 


0564 


2000 


Load TF 


0334 


31 


0521 


0562 


2002 


Load A 


0335 


34 


5000 


2100 


1516 


TF - A 


0336 


36 


2000 


0532 


2000 


(TF - A )/2 


0337 


31 


2000 


2001 


0530 


to TA 


O 5 UO 


31 


0523 


0564 


2000 


Load TF 


03^1 


31 


0521 


0562 


2002 


Load A 


03^2 


34 


3000 


2100 


1540 


TF + A 


O 5 U 5 


36 


2000 


0532 


2000 


(TF + A)/2 


03 44 


31 


2000 


2001 


0527 


to TB 


0345 


31 


0523 


0564 


0524 


TF to TE 


03 46 


30 


0534 


2100 


0547 


2 to H 


03U7 


54 


3000 


2100 


0413 


Skip to 0413 


0350 


33 


0572 


0533 


0356 


Is "SA” ^ "0” 


0351 


36 


0523 


0534 


0521 


No, TF/4 


0352 


30 


0564 


2100 


0562 


to A 


0355 


56 


0523 


0532 


0522 


TF /2 


035** 


30 


0564 


2100 


0563 


to 


0355 


34 


3000 


2100 


0565 


Skip to 0365 


0356 


31 


0523 


0564 


2000 


Load TF 


0357 


51 


0524 


0565 


2002 


Load TE 


(1 ) The test switch function would be performed by the interpolation 



selection in the overall synthesis program. 
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DIRICHLET POWER SERIES EXPANSION FOR PSEUDO-SYSTEM FUNCTION 



ADD 


INST 


M 1 


M® 


M 3 


REMARKS 


O36O 


34 


3000 


2100 


1516 


TF - TE 


0361 


56 


2000 


0534 


2000 


(TF - TE)/4 


0362 


51 


2000 


2001 


0522 


to «*. 


0363 


56 


2000 


0532 


2000 


(TF - TE)/8 


03 64 


31 


2000 


2001 


0521 


to A 


0365 


31 


0523 


0564 


2000 


Load TF 


0366 


31 


0522 


0563 


2002 


Load <*\ 


0367 


34 


5000 


2100 


1516 


TF - 


0370 


31 


0521 


0562 


2002 


Load A 


0371 


34 


5000 


2100 


1516 


TF - <*. - A 


0372 


36 


2000 


0534 


0530 


(TF - * - A )/4 


0373 


30 


2001 


2100 


0571 


to TA 


0374 


31 


0522 


0563 


2002 


Load A 


0375 


34 


5000 


2100 


1540 


TF- A 


0376 


36 


2000 


0534 


0527 


(TF - A)/4 


0377 


30 


2001 


2100 


0570 


to TB 


o4oo 


31 


0523 


0564 


2000 


Load TF 


0401 


31 


0521 


0562 


2002 


Load A 


0402 


34 


5000 


2100 


1540 


TF + A 


0405 


36 


2000 


0534 


0526 


(TF + A )/4 


0404 


30 


2001 


2100 


0567 


to TC 


0405 


31 


0522 


0563 


2002 


Load * 


0406 


34 


5000 


2100 


1?40 


TF + A + <* 


0407 


36 


2000 


0534 


0525 


(TF + A )/4 


04 10 


30 


2001 


2100 


0566 


to TD 


04ii 


31 


0523 


0564 


0524 


TF to TE 


0412 


30 


0543 


2100 


0547 


Set H to 4 


0413 


31 


0531 


0572 


2000 


Load SA 


04i4 


31 


0534 


0575 


2002 


Load "2" 


0415 


34 


5000 


2100 


1540 


SA + 2 


04l6 


31 


2000 


2001 


0531 


to SA 


0417 


31 


0535 


0576 


2004 


Load N' 


0420 


34 


5000 


2100 


1507 


Divide 


0421 


31 


2000 


2001 


0520 


Result to M 


0422 


32 


04 60 


0516 


0453 


Reset put away 


0423 


33 


0520 


2100 


0457 


If M ^ 0, 0457 


0424 


35 


2100 


2100 


0544 


No, clear J tally 


0425 


34 


0544 


2100 


0427 


If J^O, skip to 0427 


0426 


34 


5000 


2100 


0453 


No, skip to 0455 


0427 


31 


0530 


0571 


2000 


Load TA 


0430 


31 


0520 


0561 


2004 


Load M 


0431 


34 


5000 


2100 


1500 


TA X M 


0432 


31 


2000 


2001 


0530 


to TA 


0433 


31 


0600 


0601 


2000 


Load summation 


0434 


31 


0530 


0571 


2002 


Load TA 


0433 


34 


5000 


2100 


1540 


Accumulate 


0436 


30 


0455 


1505 


2006 


Build put aw ay 


0437 


32 


2006 


1677 


0442 


Build put away 
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ADD 


INST 


M 1 


M* 


M 3 


REMARKS 


0440 


55 


0442 


0532 


2006 


Build put away 


044l 


32 


2006 


1677 


0443 


Build putavay 


0442 


30 


2000 


2100 


06 l 6 


Put away 


0443 


30 


2001 


2100 


0617 


Putavay 


0444 


35 


0544 


0532 


0544 


Add 1 to J tally 


0445 


35 


0433 


0512 


0453 


Modify summation pickup 


0446 


34 


0540 


0544 


0425 


If J ^ J, repeat to 042' 


OUU 7 


36 


0547 


0532 


0547 


No, reduce H by 1 


0 U 50 


3^ 


0547 


2100 


0453 


If Hi 0, skip to 0453 


0451 


35 


0314 


0512 


0314 


No, modify h^ pickup 


0452 


34 


3000 


2100 


0314 


Repeat loop 


0453 


31 


0527 


0570 


0530 


Shift TB to TA 


0454 


31 


0526 


0567 


0527 


TC to TB 


0455 


31 


0525 


0566 


0526 


TD to TC 


0456 


34 


3000 


2100 


0413 


and repeat loop 


0457 


32 


0514 


0516 


0464 


Reset hi pickup 


0460 


31 


0600 


0601 


1310 


First term pickup 


046i 


3^ 


3000 


2100 


1011 


Convert and print H 0 (: 


0462 


30 


0532 


2100 


0544 


Set J to 1 


0463 


31 


0532 


0573 


0517 


Load FI 


0464 


31 


0620 


0621 


2000 


Pickup H, term 


0465 


32 


0544 


0532 


2106 


Extract last bit of J 


0466 


34 


2006 


2100 


0470 


If J odd, skip to 0470 


0467 


3^ 


3000 


2100 


0471 


No, skip to 0471 


0470 


32 


044l 


0537 


2001 


Change sign to minus 


0471 


31 


0517 


0560 


2004 


Load FI 


0472 


34 


3000 


2100 


1507 


Divide by F! 


0473 


31 


2000 


2001 


1310 


Load H 


0474 


34 


3000 


2100 


1011 


for conversion and pr: 


0475 


35 


0544 


0532 


0544 


Increment J by 1 


0476 


31 


0551 


0544 


2000 


Load "J" 


0477 


31 


0517 


0560 


2004 


Load FI 


0500 


34 


3000 


2100 


1500 


Multiply for (F + l)I 


0501 


31 


2000 


2001 


0517 


Putavay 


0502 


35 


0464 


0512 


0464 


Modify pickup 


0503 


3^ 


0540 


0544 


0464 


Have J terms been divid< 
appropiate factorial; 


0504 

0505 


22 

-507 


0000 0000 
Available 


0000 


Yes, halt 


0510 


00 


0000 


0001 


0000 


hj pickup reset 


0511 


35 


2000 


2005 


2000 


Putavay clear routine 


0512 


00 


0002 


0002 


0000 


Modifier for 0502 


0515 


34 


3000 


2100 


0302 


Putavay clear routine 


0514 


00 


0602 


0603 


0000 


Putavay reset 


0515 


00 


0000 


0000 


0001 


Putavay clear routine 


0516 


00 


7777 


7777 


0000 


M 1 and M® extractor 


0517 


- 0536 


Used 


for tempoary storage 




0557 


02 


0000 


0000 


0000 


Sign extractor 


0540 


00 


0000 


0000 


0000 


J 
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ADD 


INST 


M 1 


M 2 


M 3 


REMARKS 


0541 


00 


0000 


0000 


0002 


2 - for H 


0542 


00 


0000 


0000 


0003 


3 - for H 


0545 


00 


0000 


0000 


0004 


4 - for H 


0544 


00 


0000 


0000 


0000 


J tally 


0545 


00 


0000 


0000 


0000 


t storage 


0546 


00 


0000 


0000 


0001 


1 


0547 


00 


0000 


0000 


0000 


H storage 


0550 


26 


2100 


2100 


0600 


Putavay clear routine 


0551 


00 


0000 


0000 


0044 


Constant for 0476 


0552 


34 


2004 


2000 


2000 


Put away clear routine 


0553 


00 


0000 


0000 


0000 


Available 


0554 


00 


2100 


2100 


1000 


Testword for clear 


0555 ■ 


■ 0577 


Used for temporary storage 


0600 - 1000 


Used for storage of results (may be reduced to provide 



a minimum of 2 (J + l) cells. 



( 2 ) The conversion and print routine is linked in this program for 
test purposes only and should be replaced with a boot-strapping 
routine to call in next portion of the overall program. 



II 



1 
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FLOW CHART OF DIRICHLET POWER SERIES EXPANSION 






SA + "2" > SA 

SA + N' > M 



Reset PA 




i 


L 


H-l H 









\ 


v. r 

L 


mod pickup 







H >0 



Yes 




\ 


L 




TB -> TA 
TC -> TB 
TD — > TC 







B 



mi 




HLT 
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APPENDIX III 



MATRIX INVERSION 

As was noted on pages 14 and 15 of the main text of this paper, the 
final solution of Ba Hli’s method led to a system of simultaneous linear 
equations. While it is certainly possible to program the solution of a 
system of equations by the straight forward elimination technique, this 
method is awkward when many equations are involved (and for precise 
synthesis of transcendental functions, there must be many). The computer 
must do all the address computation based on N, the number of poles (or 
the degree of the denominator) and r, the relative degree between numera- 
tor and denominator. This method for n equations in n unknowns, neglect- 

2 

ing the address computations, would require n divisions to eliminate the 

2 

first unknown, (n-1) to eliminate the second, etc. In all, to reduce 
the system to the desired "one equation in one unknown 11 (we are dealing 
with square matrices), requires: 

1 i 2 . operatlons . 

i*l 

If, and only if, we have saved the intermediate equations, we may 
complete the solution in 

2 

n +n 

— 2 — operations 
or for the complete solution, 

3 9 9 

n 2 2n 

~ + n + — operations 

Such a procedure would require an inordinate amount of storage 

*The term "operation" is customarily taken to mean multiplications 
or divisions. Addition, subtractions and address modifications are not 
included in operation counts. 
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space and worse, intricate address computations. 

For simplicity of programming and best relative speed, we would 

suggest that the denominator matrix, i.e., that group of equations which 

are equated to zero, and from which we determine the q's, be solved by 

1 2 

the elimination form of the inverse, * and solve the remaining equa- 
tions for Pq, p^ P m , which incident ially will always be explicit in 

terms of the q's, by simple algebra. 

We will describe this highly efficient method of computing an in- 
verse matrix and demonstrate its simplicity and economy of operation 
by means of a simple example. Briefly, by way of introduction, we have 
the situation: 

[h[| [Qj * [h] , if we multiply by the inverse [hJ ^ 

[H] ' 1 [H] [Q] - W 1 [h] we have 

[q] « (h] ^ [hi, the desired result. 

Our immediate problem is to find a rapid means of computing the 
inverse of H. The general method of solution can be summarized: 

(I) [h l] This operation can be carried out by our 

(2) [h ^ H H ^ suggested technique operating solely on 

(3) [i H the unit matrix. 

We perform this transformation on I by bringing in H one column at 
at time, multiplying by those elements which have been transformed by the 
introduction of the columns preceding. 

^George B. Dantzig, Computational Algorithm of the Revised Simplex 
Method, USAF Project Rand Memorandum, ASTIA Document No. AD 114135, 26 
Oct., 1953. 

2 

H. M. Markowitz, The Elimination Form of the Inverse and Its 
Application to Linear Programming, Management Science, Vol. 3, No. 3, 
April, 1957. 



Ill - 2 



The actual transformation may be summarized 



(1) 


V 




j,k 


(2) 


V 






(3) 


V 


H k “ P k 




(4) 




3 *jk + P kr 


(j 


(3) 




- v 3 p kj 


(J^r) 



(j « r) 



k * 1 } 2 • « t • n 



which states (for those not versed in matrix notation), given a n x n 
matrix H (1), construct the corresponding n x n unit matrix I (2). 

Multiply the k— column of H by the unit matrix (or the partly trans- 
formed matrix) to produce the pivotal column P, , (3). Divide the ele- 

t h t h 

raents of the i — row of the unit or matrix under transformation by r~* 

element of the pivotal column, (4). Correct all other rows by subtracting 

from the element to be corrected the product of pivotal element in that 

row and the corresponding element in the row computed in Step (4), 

(5). 



A simple example will make the procedure clear, say H is 



4 2 3 
1 0 2 
2 1 2 

We start by writing down the unit matrix of corresponding order, viz., 



( 2 ) 



1 0 0 
0 10 
0 0 1 



K « 1, R a 1 
r 



( 3 ) 



1 0 0 
0 10 
0 0 1 



4" 

1 

2 



Note: This is akin to merely augment- 

ing the unit matrix by for R = 1. 
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r 

(4) 



(5) 



1, R - 1 
" 1/4 0 0 

0 10 
0 0 1 
‘1/4 0 0 

-1/4 1 0 

- 1/2 0 1 



1 u- 

j « r 



X lk + P l, 



1,2. • • • n 



X jl * X J1 ‘ P lj 1 lk (J * r * ( k “ 
Discard P. • Note: Columns to the right 

J 

of I , are not operated on, herein is the 
JK r 

great economy of this method. 



2, R 



(3) 


"1/4 


0 


0 


1/2 




-1/4 


1 


0 


-1/2 




-1/2 


0 


1 


0 


(4) 


"1/4 


0 


0 


in 




1/2 


-2 


0 


-1/2 




“1/2 


0 


1 


0 


(5) 


“ 0 


1 


0 


• 




1/2 


-2 


0 


• 




-1/2 


0 


1 


• 



V • H 2 " P 2 



J 2k 



^2k * ^2 ( seconc * element in column P) 






i - p i' 
jk 2^ 2k 



Note: Again since K e 2 <3, column 3 was 



unaffected; also note since P, 



0, I 



3k 



was unaffected. 

K ■ 3, our final transformation is: 



(3) 


0 


1 


0 


2 “ 


(4) 


"0 


1 


0 


2 ~ 


(5) 


2 


1 


-4~ 




-1/2 


-2 


0 


-5/2 




1/2 


-2 


0 


-5/2 




-4/2 -2 


5 




-1/2 


0 


1 


l/2_ 




-1 


0 


2 


l/2_ 




_-l 


0 


2_ 



and therefore 



H 



-1 



2 

-2 

-1 



1 

-2 

0 



-4 

5 

2 



III 



Check: 



[H] [H]' 1 .[l] 
^4 2 3 

1 0 2 

2 L 2 



2 L 

-2 -2 

■1 0 



-4' 

5 

2 



1 0 0 

0 L 0 

0 0 1 



The Cotal number of operations (neglecting all Incidental nulls) is: 



Multiplications to bring in each new row 

Multiplications and divisions to accomplish 
each transformation 

To solve the equations once the inverse has 
been computed 

Total product operations to complete solution 



n(n+l)(n+2)/3 

(n+l)(n)(n-l)/6 

2 



3 . 2 

n + 4n + n 



Reference to Table III-l shows the efficiency of the elimination 
form of the inverse. It is important to note that by minimizing the 
number of operations, we also minimize the propagation of round-off error. 
Von Neumann^ has determined that an approximate inverse will be found 



if 



s /4 

< .15/3 , where 2> is the base of the number system used in the 



computation, and S is the number of places. The program for the elimina- 



tion form of the inverse (appended) used a packed floating point where 
the word structure was: 

XS MMMM MMMM MQEE Legend: X « not used 

S * sign of magnitude 
M * magnitude (octal digit) 

Q ® 2 binary bits of magnitude 
& sign of exponent 
E = magnitude of exponent 



^J. Von Neumann and H. H. Goldstine, Numerical Inverting of Matrices 
of High Order, American Mathematical Society Bulletin, Vol. 53, pp. 1096. 
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This realized an S of 29, and therefore theoretically the program should 
successfully compute the Inverse for all matrices of n<23. 

The maximum round-off error to be expected is (after Von Neumann): 



.5 fi 



-(S-l) / 3 



n . „ 2 
2 * + 2n 



which for |3 ■ 2, S * 29, and n * 25 reduces to 
6 - - -000067 



Thus we can expect an accuracy of, at worst, four significant figures for 
a 25 X 25 matrix. 

If we decide to use a double precision floating point, S (for the 
CRC 102A) would increase to 66, by using the following configuration: 



XS L MMMM MMMM MMMM XS 2 MMMM MMMM MMEE 

where is sign of magnitude 
S 2 is sign of exponent 

and the round-off error now becomes (for n o 25) 



£ 



1.8125 
3.7 X 10 Ib 



« 10 



This requires a considerably greater amount of computation time 

since a floating point multiplication is now three "old" floating point 

multiplications and two floating point additions, however even for n a 60, 
_ -9 

t: < 10 . For comparative purposes, note Table III-2, The choice of 

either packed or double precision is really dependent on the order of 
the matrix expected and the requirements of the problem, it seems to this 
writer that both procedures should be available. 

It must be mentioned in passing that the method of Hotelling and 
Bodewig^ for iterating to find an improved inverse from an approximate 



A. S. Householder, Principles 
New York, 1953, pp. 80 ff. 



of Numerical Analysis, McGraw-Hill, 
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alWMMMl 



inverse and the original matrix, is readily available. However, although 

the programming of this scheme would be extremely simple, inasmuch as 

3 

the basic matrix multiplication is already present, two n multiplica- 
tions are required for each successive improvement, and furthermore, 
another n x n matrix must be stored. These requirements with regard to 
storage and computation time make the double precision floating point 

arithmetic routine more attractive, since the time to solve an n x n and 

5n 3 + 4n 2 + n 

perforin one iteration would be proportional to - X 200 for 

the packed floating point, whereas with equal storage space, the double 

n 3 + 4n 2 + n 

precision would only be proportional to x 425. 

The decision as to whether or not the simple packed floating point 
routine or the double precision was to be used would be made in the initial 
stages of the problem; i.e., when "n + r" was determined, assuming the 
number of elements was not specified (see Appendix I). 

In order to obtain a practical estimate of the accuracy of the 
matrix inversion program and the time required, the following tests were 
made : 

A symmetric 10 x 10 matrix was formed from a column of random numbers, 
the sign and position of the decimal point were also established using 
random numbers. This matrix was then inverted using the program. The 
matrix product of the original matrix and the inversion was then com- 
pared with a unit matrix. The rms . error of each element was found to 
be 431 x 10 the maximum error was 2430 x 10 

A similar test for 6x6 matrix resulted in an rms. error of 14 x 

-9 -9 

10 and a maximum error of 43 x 10 

3 

The running time was determined to be approximately .064 n minutes, 
(64 minutes for a 10 x 10, 14 minutes for a 6 x 6). 
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The flow chart and program for the elimination form of the inverse 
using simple packed floating point airthmetic is appended hereto. 
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TABLE III-l 



COMPARISON OF VARIOUS METHODS 



Method 



No. of Operations 



Word-Times (n=20)* 



S tiefel-Hestenes 


2n 


+ . • . . 


112 


X 


10 


Strict Inversion^ 


2n 3 


+ n 2 


115 


X 


io 4 


Orthogonalization^ 


2n 3 


+ n 2/ 2 


113 


X 


io 4 


2 

Elimination Form 


1/2 


. 3 . 2 . 

(n + 2n 4* n) 


31 


X 


io 4 


Back -Substitution 


1/6 


3 2 

(2n + 5n 


21 


X 


io 4 



♦Assumes product operation requires 70 word-times. 



TABLE III-2 

COMPUTATION TIME IN WORD-TIMES 

Double Precision Float 
425 
684 

♦Assuming overflow occurs 1/4 of the time. 



Operation 


Normal Float 


Packed Float 


Multiply 


80 


200 


Divide* 


112 


232 



^Householder, 

2 Von Neumann, ^ 24 ^ 
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MATRIX INVERSION PROGRAM (CRC - 102A) 



ADD INST M 1 M 2 M 3 REMARKS 



0100 


30 


0076 


2100 


0571 


0101 


30 


0077 


2100 


0366 


0102 


50 


0371 


0561 


0355 


0103 


05 


3000 


3000 


0370 


0104 


34 


3000 


2100 


2000 


0105 


35 


0377 


2100 


0106 


0106 


30 


0077 


2100 


1777 


0107 


55 


0106 


0376 


0106 


0110 


54 


0375 


0106 


0106 


0111 


35 


2100 


2100 


0375 


0112 


35 


0371 


0335 


2006 


0113 


32 


0337 


0332 


0114 


0114 


30 


0367 


2100 


0514 


0115 


35 


0114 


2006 


0114 


0116 


35 


0373 


0355 


0573 


0117 


34 


0371 


0373 


0114 


0120 


50 


0335 


2100 


0336 


0121 


35 


0334 


0366 


0365 


0122 


54 


3000 


2100 


0124 


0123 


35 


0336 


0335 


0336 


0124 


32 


2100 


0332 


0152 


0125 


35 


2100 


2100 


0373 


0126 


32 


0362 


0363 


0135 


0127 


36 


0365 


0336 


0327 


0130 


26 


2100 


2100 


0323 


0131 


34 


0327 


0326 


0133 


0132 


34 


3000 


2100 


0145 


0133 


50 


0327 


0361 


2006 


0134 


32 


2006 


0363 


01 4o 


0135 


30 


0510 


2100 


2000 


0136 


34 


2000 


2100 


0140 


0137 


34 


3000 


2100 


0145 


0140 


30 


1704 


2100 


2001 


0l4l 


54 


3400 


2100 


0400 


0142 


50 


0364 


2100 


2001 


0143 


34 


3000 


2100 


0400 


0144 


30 


2000 


2100 


0364 


0145 


55 


0325 


0335 


0323 


01 46 


35 


0135 


O36O 


0135 


0147 


35 


0327 


0355 


0327 


0150 


34 


0571 


0323 


0131 


0151 


36 


0135 


0360 


0135 


0152 


30 


0364 


2100 


0000 


0153 


35 


0373 


0335 


0375 


0154 


34 


0371 


0573 


0156 


0155 


34 


3000 


2100 


0161 


0156 


35 


0135 


O36O 


0135 


0157 


35 


0152 


0335 


0152 


0160 


34 


3000 


2100 


0127 



Store N 
Store (n + r) 

Shift N to M 1 position 
Clear cells 
0500 - 1777 
Initialize 0106 
Shift the 

contents of Channel 0 
to Channel 17 
Set RC ( row count ) to 0 
N + 1 modifier 
Initialize 01l4 
Build 
N x N 
unit 
matrix 

Set CM (column count) to 1 

1701 + (n + r) = S 

Skip 0123 initially 

Increment CM by 1 

Set put away to 0000 

Set RC to 0 

Reset pickup for ajj 

S - CM * RS( row selector) 

Clear CC (column tally) and TS 
If RS greater 1677, skip to 133 
No, skip tp 0145 
Shift RS to M 1 position 
Intalize 0140 to RS 
Pickup a^ to "A" 

If A / 0, skip to 01 40 
If A e 0, skip to 0145 
Load column element 
Transfer to multiply 
Load TS 

Transfer to add 

Store sum in TS 

Increment CC by 1 

Increment a^ * pickup by 1 

Increment Rs J by 1 

If N greater than CC, repeat 

No, compensate pickup 

Put away row-column product 

Add 1 to RC 

If N greater than RC, skip 0155 
No, skip to 0155, thence to 0l6l 
Add 1 to a^j pickup 
Increment putaway address 
Repeat loop 
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ADD 

Ol6l 

0162 

OI63 

Ol64 

OI65 

Ol66 

0167 

0170 

0171 

0172 

0173 

0174 

0175 

0176 

0177 

0200 

0201 

0202 

0203 

0204 

0205 

0206 

0207 

0210 

0211 

0212 

0213 

0214 

0215 

0216 

0217 

0220 

0221 

0222 

0223 

0224 

0225 

0226 

0227 

0230 

0231 

0232 

0235 

0234 

0255 

0236 

0237 

0240 

0241 

0242 



MATRIX INVERSION PROGRAM (CRC - 102A ) 



INST 


M 1 


M 3 


M 3 


REMARKS 


36 


0336 


0335 


0357 


CM - 1 = J 


26 


0357 


0371 


2006 


J x N 


35 


2006 


0337 


2007 


0500 + (J x N) 


35 


2100 


2100 


0523 


Clear CC 


30 


2007 


0330 


0356 


Initialize I r 


32 


0356 


0363 


0172 


in 0172 


50 


0557 


0330 


2006 


Initialize 


52 


2006 


0363 


0173 


in 0173 r 


52 


2007 


0332 


0175 


Initialize I£ in 0175 


30 


0510 


2100 


2000 


Load I r 


30 


0002 


2100 


2001 


Load P k 


54 


3600 


2100 


0400 


Ir i ? k = i; 


50 


2000 


2100 


0510 


Put away result 


55 


0323 


0335 


0323 


Increment CC by 1 


34 


0336 


0323 


0201 


If CM greater than CC, 0201 


34 


3000 


2100 


0204 


If CM = CC, 0204 


35 


0172 


0360 


0172 


Increment I r to I r 
Inc reme nt put away 


35 


0175 


0335 


0175 


34 


3000 


2100 


0172 


Continue dividing row 'r* 


35 


2100 


2100 


0325 


Set M to 0 


34 


0536 


0325 


0207 


If CM greater than M, 0207 


34 


3000 


2100 


0242 


If CM = M, 0242 


55 


0337 


0325 


0322 


Build I ' . 
Build Ijj 


30 


O322 


0350 


0324 


55 


0371 


0535 


0327 


N + 1 to TS 


30 


0335 


2100 


0373 


Initialize RC to 1 


32 


0356 


0363 


0222 


Initialize l* r pickup 


32 


2100 


0363 


0223 


Initialize P, pickup 
Initialize I pickup 

Initialize Ijj pickup 


32 


0324 


0363 


0226 


32 


0322 


0352 


0230 


34 


0375 


0336 


0222 


If RC / CM 


34 


0336 


0373 


0222 


skip to 0222 


34 


3000 


2100 


0231 


If RC * CM, skip to 0231 


30 


0510 


2100 


2000 


Load I lf 


50 


0002 


2100 


2001 


Load Pjj 


34 


3400 


2100 


0400 


Multiply 


30 


2000 


2100 


2001 


Shift for subtraction 


30 


0510 


2100 


2000 


Load 1 1 * 


34 


3200 


2100 


0400 


I’M = Mi ~ I x r x P k r 
Put away d 


30 


2000 


2100 


0510 


55 


0373 


0335 


0373 


Increment RC by 1 


34 


0327 


0373 


0236 


If N + 1 greater RC, 0236 


35 


0356 


O36O 


0356 


Increment Ii r to Is r 


35 


0325 


0535 


0325 


Increment M by 1 


34 


3000 


2100 


0205 


Repeat loop 


35 

55 


0223 

0226 


0360 

0355 


0223 

0226 


Increment P, to P^ 1 
Ixj + N - I 2 j 


35 


0230 


0371 


0230 


Increment put away 


34 


3000 


2100 


0217 


Repeat loop 


34 


0371 


0336 


0123 


Next iteration ( <x_ loop) 
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MATRIX INVERSION PROGRAM (CRC - 102A) 



ADD 


INST 


M 1 


M 2 


M 3 


REMARKS 


02^3 


35 


2100 


2100 


0336 


Set CM = 0 


0244 


35 


03U7 


2100 


0152 


Reset putavay 


0245 


30 


0367 


2100 


0000 


Put away "l" in 0000 


0246 


52 


0376 


0332 


0152 


Dummy ( see 0244 ) 


0247 


55 


035 1 * 


2100 


0155 


Mod 0155 ("out comd”) 
Transfer for N+l 8 column 


0250 


3 1 * 


3000 


2100 


0125 


0251 


35 


0346 


2100 


0152 


Reset putavay 


0252 


35 


0353 


2100 


0155 


Restore 0155 


0253 


35 


2100 


2100 


0352 


Set EC (equation count) = 0 


0254 


52 


0351 


0332 


0276 


Initialize putavay 


0255 


34 


3000 


2100 


0260 


Skip modifiers 


0256 


35 


0352 


0^35 


0352 


Increment EC by 1 


0257 


35 


0276 


0*05 


0276 


Increment putavay 


0260 


32 


0351 


0363 


0267 


Reset ho 


0261 


35 


■=<_100 


2100 


0364 


Clear TS 


0262 


30 


0352 


0361 


2006 


Shift EC 


0263 


32 


2006 


0363 


0270 


Extract into b pickup 


0261* 


5^ 


3000 


2100 


0267 


Jump modifiers 


0265 


35 


0267 


O36O 


0267 


Increment h 


0266 


56 


0270 


0360 


0270 


Decrement b 


0267 


30 


1701 


2100 


2000 


Load h 


0270 


50 


0000 


2100 


2001 


Load b 


0271 


34 


3400 


2100 


0400 


Multiply 


0272 


30 


0364 


2100 


2001 


Load TS 


0273 


54 


3000 


2100 


0400 


Add 


0274 


30 


2000 


2100 


0364 


Putavay summation 


0275 


34 


0270 


0350 


0265 


If b greater than 0, repeat 


0276 


30 


0364 


2100 


1001 


No, summation to ansver storage 


0277 


34 


0366 


0352 


0256 


If (n + r) greater EC, 0256 


0300 


55 


2100 


2100 


0345 


Clear tally, inversion completed 


0301 - 


0321 


Available for boot-strapping routines 
Denominator in Channel 0 
Numerator in Channel 1 


0322 


00 


0000 


0000 


0502 


I'ij in M 3 storage 


0323 


00 


0000 


0000 


0003 


CC storage 


032U 


00 


0502 


0000 


0000 


I ' in M 1 storage 

M storage 


0325 


00 


0000 


0000 


0003 


0326 


00 


0000 


0000 


1677 


Testword for 0131 


0327 


00 


0000 


0000 


1705 


RS storage 


0350 


00 


0000 


0000 


0030 


Shift Control 


0351 


35 


2000 


2005 


2000 


Clear Routine (2001) 


0352 


00 


0000 


0000 


3777 


M 3 extractor 


0353 


34 


3000 


2100 


0105 


Clear Routine (2003) 


033U 


00 


0000 


0000 


1701 


Constant for 01l6 


0335 


00 


0000 


0000 


0001 


Constant for Clear Routine 


0336 


00 


0000 


0000 


0000 


CM storage 


0337 


00 


0000 


0000 


0500 


Constant for 0113 etc . 
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MATRIX INVERSION PROGRAM (CRC - 102A ) 



ADD 


INST 


M 1 


M 2 


M 3 


REMARKS 


03^0 


00 


0000 


0000 


0000 


Available 


O 5 U 1 


00 


1000 


0000 


0000 


Constant for 0156 


0342 


1-0 


0100 


0000 


0000 


Not required 


05^3 


00 


0000 


0000 


0004 


Not required 


03 44 


00 


0000 


0000 


0004 


Not required 


03^5 


00 


0000 


0000 


0004 


Not required 


05 46 


30 


0364 


2100 


0000 


Reset for 0152 


0547 


36 


2100 


2000 


0000 


Reset for 0244 


0350 


00 


0000 


2100 


2001 


Testvord for 0275 


0551 


00 


1700 


0000 


1000 


Putavay reset 


0352 


00 


0000 


0000 


0001 


EC tally 


0355 


34 


3000 


2100 


0161 


Modifier for 0252 


0354 


54 


3000 


2100 


0251 


Modifier for 0247 


0355 


00 


0003 


0000 


0000 


N in M 1 


0556 


00 


0511 


0000 


0000 


hal D ly Ml 


0357 


00 


0000 


0000 


0002 


0560 


00 


0001 


0000 


0000 


Constant for 0146 


0361 


00 


0000 


0000 


0030 


Shift control word 


0362 ; 


00 


0500 


0000 


0000 


Constant for 0126 


0563 


00 


7777 


0000 


0000 


M 1 extractor 


0364 


02 


2314 


1540 


0207 


Temporary storage (TS) 


0365 


00 


0000 


0000 


1702 


S storage 


0366 


00 


0000 


0000 


0001 


(n ♦ r) storage 


0367 


00 


4000 


0000 


0001 


"1" in PackBinFltPt 


0570 


26 


2100 


2100 


0500 


Clear Routine (2000) 


0371 


00 


0000 


0000 


0003 


N storage 


0572 


34 


2004 


2000 


2000 


Clear Routine (2002) 


0575 


00 


0000 


0000 


0003 


RC tally 


0374 


00 


2100 


2100 


l4oo 


Clear Routine (2004) 


0375 


00 


0077 


2100 


1777 


Testvord for 0105 


0376 


00 


0001 


0000 


0001 


Constant for 0104 


0577 


30 


0000 


2100 


1700 


Not required 


Channel 4 contains packed floating point 


arithmetic routines 


o4oo 


30 


3000 


0445 


2006 


Prepare 


0401 


32 


2006 


0475 


0455 


exit 


0402 


30 


2001 


2100 


2007 


Store 2 n< * operand 


0403 


52 


2000 


0454 


2102 


Extract sign of exp 1 


0404 


27 


2002 


0451 


2002 


Shift to operating location 


0405 


32 


2000 


0453 


2002 


Extract mag of exp 1 


0406 


32 


2000 


0444 


2105 


Extract sign and mag of mag 1 


0407 


31 


2002 


2003 


2000 


Shift to 2000,2001 


04l0 


32 


2007 


0434 


2102 


Extract sign of exp 2 


04ll 


27 


2002 


0451 


2002 


Shift to operating location 


0412 


32 


2007 


0436 


2002 


Extract mag of exp 2 


0415 


32 


2007 


0444 


2103 


Extract sign and mag of mag 2 


04l4 


36 


2006 


0455 
III - 


2006 

1-4 


Pickup entry command 



MATRIX INVERSION PROGRAM (CRC - 102A) 

ADD INST M 1 M 2 M 3 REMARKS 



0415 


30 


2006 


0437 


2006 


04i6 


32 


2006 


046l 


0417 


0417 


32 


0273 


046l 


2107 


0420 


54 


2007 


0400 


0440 


0421 


33 


2002 


2000 


0467 


0422 


36 


2002 


2000 


2006 


0425 


30 


2003 


2006 


2007 


0424 


35 


2001 


2007 


2001 


0425 


37 


2001 


3000 


0462 


0426 


31 


2000 


2001 


2002 


0427 


27 


2002 


0452 


2104 


0430 


32 


2002 


0436 


2004 


0431 


32 


2004 


0 476 


2003 


0432 


30 


2003 


2100 


2000 


0433 


34 


3000 


2100 


0274 


0434 


00 


0000 


0000 


0100 


0435 


00 


0000 


0000 


0001 


0436 


00 


0000 


0000 


0077 


0437 


00 


0000 


0000 


0030 


o44o 


34 


2007 


0450 


0445 


044i 


36 


2100 


2003 


2003 


0442 


34 


3000 


2100 


0421 


0443 


02 


0000 


0000 


0030 


0444 


02 


7777 


7777 


7600 


0445 


34 


2007 


0474 


0454 


0446 


55 


2000 


2002 


2000 


0447 


25 


2001 


2005 


2001 


0450 


34 


3200 


2100 


0426 


0451 


00 


0000 


0000 


0037 


0452 


02 


0000 


0000 


0037 


0453 


00 


0000 


0000 


0011 


0454 


36 


2000 


2002 


2000 


0455 


23 


2001 


2003 


2001 


0456 


37 


2001 


3000 


0462 


0457 


54 


3000 


2100 


0426 


0460 


02 


0000 


0000 


0001 


046i 


00 


3777 


0000 


0000 


0462 


27 


2001 


0460 


2001 


0463 


30 


2001 


0460 


2001 


0464 


27 


2001 


0477 


2001 


0465 


33 


2000 


0477 


2000 


0466 


34 


3000 


2100 


0426 


0467 


56 


2000 


2002 


2006 


0470 


30 


2001 


2006 


2007 


0471 


30 


2002 


2100 


2000 


0472 


35 


2005 


2007 


2001 


0473 


37 


2001 


3000 


0462 


0474 


34 


3400 


2100 


0426 


0475 


00 


0000 


0000 


7777 


0476 


00 


0000 


0000 


0177 


0477 


00 


0000 


0000 


0001 



Shift to M 1 
Set OU17 
Extract control 
If control ^ 3000, 0440 
No, operations 
for normal 
floating point 
addition 
routine 

Store in 2002, 2003 
Shift sign of exp for pack 
Extract mag of exp for pack 
Pack exponent 
Packed result to 2000 
Exit to main program 
Extractor for sign of exp 
Constant for 04l4 
Extractor for mag of exp 
Shift control word 
If control ^ 3200 0445 

Change sign of 2 n ” operand 
Transfer to addition 
Shift control word 
Extractor for magnitude 
If control ^ 3400, 0454 
No, operation for normal 
floating point 
multiplication 
Shift to pack commands 
Shifter for sign of exp 
Extractor for mag of exp 
Operations for 
normal floating point 
division 

Transfer to pack commands 
Constant for 0462 
M 1 extracotor 
Normal 
correct 
overflow 
procedure 

Transfer to pack commands 
Operations for 
normal floating 
point addition 
when 2002 > 2000 
see 0421 

Transfer to pack commands 
M 3 extractor 

Sign and exponent extractor 
Constant 
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FLOW CHART OF MATRIX INVERSION AND SOLUTION 



C ''Matrix" [H] ) 



1 



Clear 0500-1777 
p}— 1700 
Set RC - 0 
Set N+l in modifier 
Initialize "1" PA 



*<D 



This N is the 
order of the 
matrix, viz. n 



Add N+l to 
"1" PA 
Cooxaand 




■>Add 1 to CM 



11 1 H 



ij 



(i-j) 

add 1 to RC 



Yes 




Set 1 +CH 
1701 + n + r-^S 



Unit Matrix 



n + r 



3 




Reset PA to 0000 
Set RC to 0 

Reset a^ to a^ pickup 



S - CM-+RS 
Clear CC 
Clear TS 



Setup to 
multiply 
in rows 



Row selection 







r 



No 



^ 








r 











add 1 to EC 
Mod PA 

X 



1 

Reset h E( , 

Clear TS 
Shift EC 

Initialize pickup 



APPENDIX IV 



LAG RANG IAN INTERPOLATION METHOD FOR 
THE SOLUTION OF ALGEBRAIC EQUATIONS 1 

The Lagraagian interpolation method described by D. E. Muller was 
selected for the solution of large degree polynominals . This procedure 
was programmed, with certain modifications, and tested on the CRC-102A 
computer. By removing the roots, as found, by synthetic division, rather 
than by the algorithm recommended by Muller, accuracies of at least eight 
decimal places were obtained for all roots of equations up to degree 21. 
This was considerably better than the accuracies reported by Muller, and 
incidentally, probably entailed considerably more iteration. 

-13 

A graph of solution time vs. degree of equation, for £ * 10 , is 

presented in Figure IV-1. Convergence to a root generally requires be- 

-13 

tween 7 and 14 iterations (for the £ of 10 ). The iterative loop is 

completed in about 1.5 minutes, varying slightly (not more than 10%) 

according to the degree of the equation. In making these tests, it was 

found that when equations had roots of the same modulus, there is a pro- 

24 

nounced slow down in the rate of convergence. The equations x +1*0, 

20 

and x +1*0 required about 12.5 minutes per root, whereas the more 
general polynominal, (equation A under Results of Tests) only required 
about 9.4 minutes per root. Round-off error was not present in the 

2 0 24 

x +1 equation and affected only two roots in the solution of x +1*0 
within the eighth decimal place. 

^D. E. Muller, Solving Algebraic Equation by an Automatic Computer, 
Mathematical Tables and Other Aids to Computation, Vol. X, No. 56, Oct., 
1956, pp. 208-215. 
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Thus far the program has failed only once. When solving an equation 
with a real double negative root and a pair of complex roots of about 
the same modulus , viz., 4; the rate of convergence was normal until about 
four decimal places had been established and decreased so as to become 
barely perceptible. The program did, however, remove all roots interior 
to modulo 4 without difficulty. 

The basic outline of the Lagrangian interpolation method is as 
follows : 

Successive iterations toward any particular root are obtained by 
finding the nearer root of a quadratic whose curve passes through the 

points ; z^, f(z^)j ^ z i-l^* Z i~2 9 ^ z i-2^* 

The Lagrangian formula may be written as £f(z)J * b ^z ^ + b^z + b^; 

where the b's satisfy the system of equations: 

Vi 2 + Vi + b 2 “ f( V 
b 0 Z i-l + Vi-1 + b 2 “ f(z iV 
Vi-2 + Vi-2 + b 2 E f(z i-2> 

and these equations, of course, become almost one in the same as the 

iterative process converges, depending on the accuracy which we specify. 

We are to construct a polynominal (quadratic) whose curve passes 

2 

through the three points. By the Lagrangian interpolation formula , we 
write : 



r i < z - ” z i -2 ) 

L. f(z)L f( 2j ) t ■ v- — — r 

1 L J 1 < z i ' z i-l )(z i " Z i-2> 



+ f ( z <_,) 



( Z - ^Hy^) 
i-i' (« 1 . l - « 1 )(« 1 . 1 - y 2 ) 






(z - z t ) ( z - z i l ) 
l * 2 ' (l i-2 • z i )<z l-2 • z l-l> 



W. E. Milne, Numerical Calculus, Princeton University Press, Princeton, 
N. J., 1949, for example. 
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If we make the following simplifying substitutions: 



z - » h 



z i • z i-i * h i 



h/h^ * X 

h i /h i-i 



'i-l 



- z 



i-2 



i-l 




A 



i 



we may express the coefficient of f(z^) as: 

% Ai 1 , + + i 

<z ■ z^Kz - z 1-2 ) _ + n < z - y 2 > 

(Z i ' Z l-l )(Z t • Z i-2 ) " L <Z i ' Z l-2 ) 



since 



1 z i-l " z i-2 



\ z i ' Vi 



, and 



i + X + 



z - z 



i-2 



z i • vi \ 



^ * X 






fX x 



a 'A xr A 



z - z 



i-2 



Z i-1 " Z i-2 



z - z 



i-2 



Z i " Z i-2 



and finally 



; therefore. 



(z - z. )(z - z, . ) 



(1 +X) (1 +XX.A' 1 ) - 7 — T7 — — 7 the coefficient of f (z ), 

11 U i z i-l MZ i ‘ Z i-2 ' 1 



which when expanded is: 



^ X i Ai 1 + X\i Aj^ 1 + X + 1 = ^ \ A^ 1 + \&\\ +Ap + i (QED). 



By similar manipulation, the coefficients of f(z^_^) and f(z^_^) cart 
found and the desired quadratic in \ is: 
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A single iterative step is obtained by letting z ^ be such a value 
of z that [f(z)] * 0. We may solve the quadratic inversely to obtain 



the iterative process* 

All operations are in complex numbers and in floating point arithmetic. 
A root is taken when 



For proof of convergence of the process the reader is referred to Dr. 
Muller*s paper* 




from which we obtain the new z^ and repeat 
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RESULTS OF TESTS 



Equation A 



Z 16 - 9.0 Z 15 + 34.75 z 


14 - 102 


13 12 11 

.25 z + 264.0 z - 479.0 z 


10 9 

+ 771.5 z - 1154.5 z 


+ 978.0 


z 8 - 1225.0 z 7 + 912.75 z 6 


+ 254.75 z 5 + 1411.0 z 4 


+ 1319. 


0 z 3 + 1022.5 z 2 + 470.5 z 


-44 

+ 105.0 £ a 2 (. 10 


- 13 > 




Total solution time: 2 


hours , 


34 minutes, 10 seconds on CRC-102A. 


Roots : 






-0.284460054 


± j 


0.598759363 


-0.376766336 


± j 


0.188840569 


-0.101279181 


± J 


1.429775042 


-0.044975589 


± j 


1.738162299 


+1.999999996 






-0.011258189 


± j 


1.933289317 


+2.500000009 






-0.181260648 


± j 


1.034830851 


+3.000000001 






+3.499999991 
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RESULTS OF TEST (Continued) 



Equation B 



30 

x 



+ 1 



0 



6*2 



-44 



(- 10 



-13 



) 



Total solution time: 6 hours, 23 minutes, 30 seconds on CRC-102A 

Roots : 



- 0.000000000 


± J 


0.999999999 


-0.951056516 


+ J 


0.309016994 


+0.866025403 


± J 


0.500000000 


-0.587785252 


± J 


0.809016994 


+0.994521895 


1 J 


0.104528463 


-0.207911690 


± j 


0.978147601 


-0.994521895 


± J 


0.104528463 


+0.406736637 


± J 


0.913545463 


+0.207911688 


+ i 


0.978147598 


+0.743144829 


± J 


0.669130603 


-0.743144826 


± J 


0.669130607 


+0.587785258 


± J 


0.809017000 


-0.866025404 


+ j 


0.499999998 


+0.951056517 


± j 


0.309016993 


-0.406736643 


± J 


0.913545458 
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LAGRANGIAN INTERPOLATION METHOD FOR HIGH DEGREE POLYNOMIALS 



ADD INST M 1 M 2 M 3 REMARKS 



Channels 0 and 1 are used temporary storage. 

Channel 2 contains a decimal to binary floating point conversion 
which while not needed in overall program is Included here. 



0200 


26 


1540 


0250 


2005 


0201 


55 


0246 


2005 


0255 


0202 


55 


0251 


2100 


0210 


0205 


55 


0255 


2100 


0211 


0204 


52 


2100 


0252 


0270 


0205 


05 


5000 


5000 


5000 


0206 


50 


0254 


2100 


0265 


0207 


52 


0254 


0255 


2004 


0210 


50 


0016 


2100 


0250 


0211 


50 


0017 


2100 


0247 


0212 


54 


0250 


2100 


0214 


0215 


54 


5000 


2100 


0251 


0214 


52 


0250 


0255 


0227 


0215 


52 


0250 


0265 


2005 


0216 


26 


2004 


2005 


2006 


0217 


55 


2006 


2001 


2001 


0220 


26 


2004 


0257 


2004 


0221 


50 


0250 


0262 


0250 


0222 


54 


0250 


2100 


0215 


0225 


51 


2000 


2001 


2000 


0224 


55 


2000 


0261 


2000 


0225 


54 


0247 


2100 


0255 


0226 


54 


5000 


2100 


0267 


0227 


00 


0000 


0000 


0000 


0250 


00 


0000 


0000 


0000 


0251 


54 


0247 


2100 


0254 


0252 


54 


5000 


2100 


0267 


0255 


00 


0016 


2100 


0250 


0254 


52 


0247 


0255 


0227 


0255 


52 


0247 


0265 


2005 


0256 


25 


2005 


0265 


2005 


0257 


50 


0265 


0256 


0265 


0240 


52 


0247 


0265 


2005 


0241 


25 


2005 


0265 


2005 


0242 


54 


0260 


0265 


0257 


0245 


54 


2001 


2100 


0266 


0244 


51 


2002 


2005 


2000 


0245 


54 


5000 


2100 


0267 


0246 


00 


0002 


2100 


0250 


0247 


00 


0000 


0000 


0000 


0250 


00 


0002 


0000 


0000 


0251 


50 


0000 


2100 


0250 


0252 


00 


0000 


0000 


7777 



Build and 

store testword 
Reset pickup 
Reset pickup 
Reset put away 
Clear buffer 
Reset extractor 
Set 2004 to 1 
Pickup integer 
Pickup fraction 
Is integer ^ 0? 

No, test fraction 
Yes, extract sign 
Extract least sig. digit 
L.S.D. x (2004) 

Accumulate products 
Increase 2004 by power of 10 
Shift decimal word 4 right 
If conversion not complete, 0215 
If complete normalize 
and correct exponent 
Is fraction ^ 0? 

No, skip to put away 
Sign storage 
Integer storage 
I * 0, Is fraction ^ 0? 

If both are 0, put away 0 
Testword 

Extract sign when fraction only 

Extract least sig. frac. dec. dig. 

Divide by 10 

Shift extractor 

Extract next digit 

Divide by 10 

If conversion not complete, 0257 
If integer 0, transfer to add 
If integer was 0, normalize 
and put away 
Testword addend 
Fraction storage 
Testword build constant 
Integer pickup reset 
M 3 extractor 
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LAGRANGIAN INTERPOLATION METHOD FOR HIGH DEGREE POLYNOMIALS 
ADD INST M 1 M 2 M 3 REMARKS 



0255 


50 


0001 


2100 


02^7 


025 k 


00 


0000 


0000 


0017 


0255 


02 


0000 


0000 


0001 


0256 


00 


0000 


0000 


000 k 


0257 


00 


0000 


0000 


0012 


0260 


00 


7k00 


0000 


0000 


0261 


00 


0000 


0000 


OOkk 


0262 


02 


0000 


0000 


000k 


0263 


00 


5000 


0000 


0000 


026 k 


00 


0000 


0000 


0002 


0265 


00 


7k00 


0000 


0000 


0266 


3^ 


3000 


2100 


15 kO 


0267 


32 


0227 


0255 


2001 


0270 


30 


2000 


2100 


00l6 


0271 


36 


0270 


0255 


2007 


0272 


32 


2007 


0252 


0273 


0273 


30 


2001 


2100 


0015 


027 k 


35 


0210 


0250 


0210 


0275 


55 


0211 


0250 


0211 


0276 


35 


0270 


026 k 


0270 


0277 


3k 


0233 


0210 


0205 



Fraction pickup reset 
Extractor reset 
Sign extractor and -1 
Shift control - 1* left 
"Power of 10" 

Extractor testvord 
To correct exponent 
Shift control ~ 4 right 
10 for division 
2 

Extractor working storage 
Transfer to addition 
Extract sign into mag portion 
Putavay exponent 
Build magnitude 
putavay 

Putavay magnitude 
Modify integer pickup 
Modify fraction pickup 
Modify putavay 

Have all numbers been converted? 



Channels 3 through 7, and 1000 through 1012 contain main program 
for Lagrangian Interpolation Method; in the program detailed 
here Channels 10 and 11 contain a binary floating point to 
decimal conversion routine which would not be required for 
the overall program. Channel 13 is used for tempoary storage 
and Channels Ik through l6 contain the floating point arith- 
metic routines for complex numbers. 
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LAGRANGIAN INTERPOLATION METHOD FOP HIGH DEGREE POLYNOMIALS 
ADD INST M 1 M 2 M 3 REMARKS 



0300 


30 


1340 


2100 


1372 


0301 


05 


3000 


3000 


0000 


0502 


04 


3000 


3000 


0100 


0503 


35 


05 Cl 


1375 


0501 


0304 


35 


0302 


1376 


0302 


0305 


34 


1374 


0502 


0301 


0506 


36 


0301 


1437 


0301 


0307 


36 


0302 


1437 


0302 


0310 


30 


1531 


2100 


2000 


0311 


36 


2100 


1537 


2001 


0312 


31 


2000 


2001 


1300 


O 513 


51 


1333 


1357 


1320 


0314 


31 


1331 


1357 


1301 


0315 


31 


1333 


1337 


1321 


0316 


31 


1333 


1337 


1302 


0317 


51 


1535 


1337 


1322 


0320 


35 


2100 


2100 


1303 


0321 


36 


2100 


1337 


1344 


0322 


31 


1333 


1337 


1325 


0323 


26 


1372 


1466 


2007 


0324 


26 


2100 


2100 


1434 


0325 


26 


2100 


2100 


1435 


0526 


35 


0640 


2007 


0327 


0327 


31 


010-2 


0105 


2000 


0330 


31 


2000 


2001 


1307 


0331 


31 


1533 


1357 


1327 


0332 


56 


0327 


1473 


2007 


0333 


32 


2007 


1464 


0554 


0534 


31 


0070 


0071 


2002 


0335 


54 


3000 


2100 


1540 


0336 


56 


0334 


1473 


0534 


0337 


54 


0640 


0334 


0541 


0340 


34 


3000 


2100 


0334 


0341 


34 


l4 15 


2100 


0547 


0342 


31 


2000 


2001 


1434 


0343 


36 


0527 


1466 


0344 


0544 


31 


0100 


0101 


2000 


0345 


36 


0344 


1475 


2007 


0546 


34 


3000 


2100 


0333 


0347 


31 


2000 


2001 


1435 


0350 


31 


1434 


l4?5 


2002 


0351 


34 


3000 


2100 


1540 


0352 


31 


2000 


2001 


1306 


0353 


51 


1335 


1557 


1526 


0354 


31 


1434 


1475 


2000 


0355 


31 


1435 


1476 


2002 


0356 


34 


3000 


2100 


1516 


0557 


31 


2000 


2001 


1305 



N to N' 

Transfer 
the contents 
of Channel 0 
to 

Channel 1 
and 
reset . 

Initial setup 
for z i= 2 

Initial setup 
for 

Initial setup 
for zj 

Initial setup 
for X; 

- 1/2 + JO 

Commence pickup build 
Clear 1434, 1475 
Clear 1 * 435 , U 76 
Set pickup for A 
Pickup A 0 
Ao is tWof f (zi ) 

<JWf (zi ) * 0 
Build A n pickup 
and enter in 0354 
A n pickup command 
Transfer to floating addition 
Modify A n pickup for A n-g 
If sum complete, putavay 
No, repeat to 0334 
Has even sura putavay been made? 
No, even sum to 1434, 1475 
Build A^^ pickup 
A n-1 pickup 

Set An pickup for odd terms 
and repeat back for odd sum 
Putavay odd sum in 1435, 1476 
Load even sum 
Odd + even for 

f ( z i -1 ) 

and clear imaginary 
Load even sum 
and odd sum 
for subtraction 
to form 



IV 
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LAGRANGXAH INTERPOLATION METHOD FOR HIGH DEGREE POLYNOMIALS 



ADD INST M 1 M 2 M 3 REMARKS 



0360 


51 


1533 


1357 


1325 


0361 


31 


1303 


13bb 


2000 


0362 


31 


1323 


136b 


2002 


0363 


31 


1331 


1335 


2004 


0364 


31 


1535 


1337 


2006 


0365 


31* 


3000 


2100 


1410 


0366 


31 


2000 


2001 


1304 


0367 


31 


2002 


2003 


1324 


0370 


31 


1306 


13^7 


2004 


0371 


31 


1326 


1367 


2006 


0372 


34 


3000 


2100 


1607 


0373 


01* 


3000 


3000 


1330 


0374 


31 


1305 


1346 


2000 


0575 


31 


1325 


1566 


2002 


0376 


31 


1303 


13l*l* 


2004 


0377 


31 


1323 


136 b 


2006 


o4oo 


34 


3000 


2100 


1607 


o4oi 


01* 


3000 


3000 


1510 


0402 


31 


1507 


1350 


2004 


0403 


31 


1527 


1570 


2006 


0404 


34 


3000 


2100 


1410 


0^05 


31 


1330 


1371 


2004 


0406 


31 


1332 


1373 


2006 


0407 


31* 


5000 


2100 


l4l6 


04l0 


31 


1301* 


13l*5 


2004 


o4n 


31 


1321* 


1365 


2006 


0412 


3b 


5000 


2100 


1607 


04l3 


31 


1303 


13l*l* 


2004 


o4i4 


31 


1323 


1361* 


2006 


0415 


3b 


5000 


2100 


1607 


0416 


31 


1307 


1350 


2004 


0417 


31 


1527 


1370 


2006 


0420 


35 


2001* 


0261* 


2004 


0421 


35 


2006 


0264 


2006 


0422 


3b 


3000 


2100 


1607 


0423 


01* 


3000 


3000 


1440 


0424 


05 


3000 


3000 


1330 


0L25 


31 


1301* 


13l*5 


2004 


0426 


51 


1321* 


1365 


2006 


0427 


31* 


3000 


2100 


1607 


OU3O 


01* 


3000 


3000 


1330 


0U31 


05 


3000 


3000 


1310 


01*32 


31 


1303 


131*4 


2004 


OU33 


31 


1323 


1364 


2006 


0L3I* 


3b 


3000 


2100 


1607 


01*55 


ou 


5000 


3000 


1310 


01*36 


31 


1303 


13l*l* 


2000 


01*37 


31 


1323 


1364 


2002 



and clear imaginary 
Load 

Xi 

Load 
1 + JO 

Transfer to cpx. add. 

Put away 

Load 

f(z 1= l) • 

Transfer to cpx. mpy . 

Store in B 
Load 

Load 

>i 

Transfer to cpx. mpy. 

Store in A 
Load 

**(* i) 

f ^al “ 2 ^i + f ^ z i ) 

_ f ( z i~i) A i 

Lf ( z i=2M i “ + 

Load 

fttL A. 

Load 

mV A 

Load 
f fz i ) 

and "multiply 
by 4 

4f X j f (z) » 4*AC 

Store in C 
Load f )• A 
Load 

A i , . a 
f fz 1= l) A i 
Store In B 
f ( z i-2 ) X i f rom A 
Load X i 

Multiply for f (Zi„2 ) X* 2 

Store in A 

Load 

Xi 
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LAGRANGIAN INTERPOLATION METHOD FOR HIGH DEGREE POLYNOMIALS 



ADD 


INST 


M 1 


M 2 


M 3 


REMARKS 


0440 


31 


1304 


1345 


2004 


Load 


o44i 


31 


1323 


1365 


2006 


^ i 


0442 


34 


3000 


2100 


1410 


( A 1 ) 


0443 


31 


1307 


1350 


2004 


Load 


0444 


31 


1327 


1370 


2006 


ffzi) 


0445 


34 


3000 


2100 


1607 


f(z*) (>, + 


0 446 


31 


1310 


1351 


2004 


Load 1 1 


0447 


31 


1312 


1353 


2006 




0450 


34 


3000 


2100 


1410 


f (zi ) T A, + Aj ) + f (z 


0451 


31 


1330 


1371 


2004 


Load 


0452 


31 


1332 


1373 


2006 


f(zj ,) 


0*53 


34 


3000 


2100 


l4l6 


f(Zi-2^Ai “ 


0454 


04 


3000 


5000 


1330 


Result =• g^; store in 


0455 


31 


2000 


2001 


2004 


Load B for complex 


0^56 


31 


2002 


2003 


2006 


multiplication, i.e. 


0457 


34 


3000 


2100 


1607 


gi 2 OR B 2 


0460 


31 


l44o 


1401 


2004 


Load 4AC 


046i 


31 


1442 


1403 


2006 


B 2 - 4AC 


0462 


34 


3000 


2100 


l4l6 


Store in C 


0463 


04 


3000 


3000 


l44o 


Load 


0 464 


31 


2000 


2001 


2004 


B 2 - 4AC 


0465 


30 


2002 


2100 


2006 


for conjugate 


0 466 


36 


2100 


2003 


2007 


multupllcation 


0467 


34 


3000 


2100 


1607 


(b 2 -4ac Hb 2 - 4ac r = 


0U70 


34 


3000 


2100 


1551 


Compute r 


0471 


31 


2000 


2001 


1434 


Store r in TS X 


0472 


34 


3000 


2100 


1551 


Compute -fr" 


0473 


31 


1434 


1475 


2004 


Load r 


0474 


31 


2000 


2001 


1434 


Store -fir in TS. 


0475 


31 


1440 


1401 


2000 


(B 2 = 4AC) = x 


0476 


34 


3000 


2100 


1507 


x/r * cos 0 


0477 


31 


2000 


2001 


1435 


Store cos & in TS 2 


0500 


31 


1331 


1337 


2002 


Load 1 


0501 


34 


3000 


2100 


1540 


1 + coe 0 


0502 


36 


2000 


1331 


2000 


(1 + cos 9 )/2 


0503 


34 


3000 


2100 


1551 


^(1 + cos Q )J2 * cos < 


0504 


31 


1434 


1475 


2004 


Load -fr 


0505 


34 


3000 


2100 


1500 


-Jr" (cos 9/2) = ^ 


0506 


31 


2000 


2001 


1436 


<*’ to TS 3 


0507 


31 


1435 


1476 


2002 


Load cos S 


0510 


31 


1331 


1337 


2000 


Load 1 


0511 


34 


3000 


2100 


1516 


1 - cos 0 


0512 


36 


2000 


1331 


2000 


(1 - cos0)/2 


0513 


34 


3000 


2100 


1551 


-((l - cos 9 )/2 = sin t 


0514 


31 


1434 


1475 


2004 


Load -fr 


0515 


34 


3000 


2100 


1500 


-fr (sin 0/2) * a' 


0516 


31 


2000 


2001 


1434 


£>’ to TSi r 


0517 


30 


2000 


2100 


1451 


Putavay (3 0 


0520 


30 


2000 


2100 


1453 


and exps 



IV - 1 - 5 



LAGRANOIAN INTERPOLATION METHOD 



ADD 


INST 


M 1 


M® 


M 3 


0521 


30 


11*36 


2100 


ll*50 


0522 


30 


11*36 


2100 


ll*52 


0523 


32 


l**75 


1575 


2107 


052L 


50 


2007 


2100 


ll*6l 


0525 


36 


2100 


2007 


11*63 


0526 


32 


1*77 


1575 


2107 


0527 


33 


ll*03 


2100 


0553 


0530 


30 


2007 


2100 


11*62 


0531 


36 


2100 


2007 


ll*60 


0532 


3* 


3000 


2100 


0535 


05 33 


30 


200? 


2100 


1*60 


05 3 1 * 


36 


2100 


2007 


ll*62 


0535 


05 


3000 


5000 


1330 


0536 


31 


ll*50 


1*60 


2001* 


0537 


31 


ll*51 


l**6i 


2006 


05*0 


3* 


5000 


2100 


l*l6 


05*1 


31 


2000 


2001 


2004 


05*2 


30 


2002 


2100 


2006 


05*3 


36 


2100 


2003 


2007 


05** 


3* 


5000 


2100 


1607 


05*5 


31 


2000 


2001 


1*3* 


05*6 


05 


3000 


3000 


1330 


05*7 


31 


ll*52 


11*62 


2001* 


0550 


51 


1*53 


1*63 


2006 


0551 


3* 


3000 


2100 


1*16 


0552 


31 


2000 


2001 


2001* 


0553 


30 


2002 


2100 


2006 


055* 


36 


2100 


2003 


2007 


0555 


5* 


3000 


2100 


1607 


0556 


51 


2000 


2001 


11*36 


- 0557 


05 


3000 


3000 


1330 


0560 


33 


ll*3l* 


11*36 


0566 


0561 


33 


11*36 


ll*3l* 


0563 


0562 


3* 


1*75 


1*77 


0566 


0563 


31 


11*51 


l*6l 


2006 


056* 


31 


ll*50 


1*60 


2001* 


0565 


3>* 


3000 


2100 


0570 


0566 


31 


11*53 


IL63 


2006 


0567 


31 


1*52 


]1*62 


2001* 


0570 


3* 


3000 


2100 


1*10 


0571 


0* 


3000 


5000 


1530 


c 572 


31 


1307 


1350 


2000 


0573 


31 


1527 


1370 


2002 


0574 


35 


13d* 


1531 


200* 


0575 


36 


2100 


13*5 


2005 


0576 


35 


1321* 


1331 


2006 


0577 


36 


2100 


1365 


2007 



IV - 1 - 6 



FOR HIGH DEGREE POLYNOMIALS 
REMARKS 



Put away c< 0 
and exp s 

Set /3 0 mag positive 
and 

_/3 x mag negative 
Clear mag sign bit 
If «fU,(B 2 - 4AC) * 0, 0533 
No, putavay °<i positive 
and <>< 0 negative 
Skip to 0535 
Putavay «>< 0 positive 
and =*! negative 
Load g i 
Load ®< 0 
Load (3 q 
S ubtract 
Load for 
conjugate 
multiplication 
Multiply and 
store in TSj 
Load gj^ 

Load 

Load 

Subtract 
Load for 
conjugate 
multiplication 
Multiply and 
store in TS 3 
Load g. 

If TSi ^ TS 3 , 0566 
If TS 3 * TS X , 0563 
If TSj. ^TS 3 , 0566 
Load 
and 
or 

load 
a nd &■> 

B +4 V ? -LAC 
Store in B 
Load 
f(z, ) 

Load 



-2&J f(z t ) 



LAGRANGIAN INTERPOLATION METHOD FOR HIGH DEGREE POLYNOMIALS 



ADD 


INST 


M 1 


M 2 


M 3 


REMARKS 




0600 


34 


3000 


2100 


1607 


Multiply 




0601 


31 


1330 


1371 


2004 


Load 




0602 


31 


1532 


1373 


2006 


B + -s/B 2 -4 AC 




0603 


34 


3000 


2100 


1650 


Divide for,\. . 




06o4 


04 


5000 


3000 


1330 


Store in B 




0605 


31 


1302 


15^3 


2000 


Load 




0606 


31 


1322 


1363 


2002 


Loai 




0607 


31 


1301 


1342 


2004 




0610 


31 


1321 


1362 


2006 


z i~l 




06ll 


34 


3000 


2100 


l 4 i 6 


z i " z i-l ■ h i 




0612 


51 


1330 


1371 


2004 


Load 




0613 


31 


1332 


1373 


2006 


Loai +1 ^ h ^ = hi+1 




06l 4 


34 


JOOO 


2100 


1607 




0615 


51 


1302 


13^3 


2004 




06l6 


31 


1322 


1363 


2006 


Z i 




0617 


34 


3000 


2100 


l 4 lo 


z i + h i+l * z i+l 




0620 


Ob 


3000 


3000 


1310 


Store in A 




0621 


31 


1302 


13^3 


2004 


Load 




0622 


31 


1522 


1363 


2006 


z i 




0623 


34 


3000 


2100 


l 4 i 6 


z i+l “ z i ■ 




0624 


31 


1302 


13^3 


2004 


f\7 z i * s 




0625 

0626 
0627 


31 

34 

33 


1522 

5000 

14^4 


1363 

2100 

2000 


2006 

1650 

0631 


If 6 ^ 1 for real part, 


0631 


0650 


34 


3000 


2100 


0632 


No, evaluate 




0631 


33 


l 4^4 


2002 


0725 


If 6 - 0 for imag part, 


0725 


0652 


26 


1372 


1466 


2007 


Build and 




0633 


35 


1470 


2007 


0656 


store testvord 




06 3L 


30 


1470 


2100 


2007 


Reset A^ pickup 
in 0646 




0635 


35 


2007 


0264 


0646 




O636 


31 


1310 


1351 


2004 


Load 




0637 


31 


1312 


1353 


2006 


z i+l 




o 64 o 


31 


0100 


0101 


2000 


Load 




o 64 i 


51 


1333 


1337 


2002 


An 




0642 


52 


0637 


1636 


1635 


Set exit in cpx. rapy. routine 


0643 


34 


3000 


2100 


l6ll 


Transfer into cpx. mpy. 




064 U 


31 


2002 


2003 


1602 


Temp, store imag part 




0645 


32 


0677 


1636 


1550 


Set exit in add routine 




0646 


31 


0104 


0105 


2002 


Load A. 




0647 


34 


3000 


2100 


1542 


Transfer Into add routine 


0650 


31 


1602 


1643 


2002 


Load imag part for cpx. 


mpy. 


0651 


31 


1310 


1351 


2004 


Reload z 1+ ^ 




0652 


31 


1312 


1333 


2006 


for next term 




0653 


35 


0646 


1466 


0646 


Modify A^ pickup 




0654 


34 


0646 


0656 


0660 


If A 0 not taken, 




0655 


34 


3000 


2100 


l6ll 


repeat back 




0656 


31 


0104 


0105 


2000 


Testvord storage 




0657 


00 


0000 


0000 


0644 


Constant for 0&*2 




0660 


04 


3000 


3000 


l 44 o 


Store f(Zi+i) 





IV - 1 - 7 



LACRANGIAN INTERPOLATION METHOD FOR HIGH DEGREE POLYNOMIALS 



ADD INST M 1 M 2 M 3 REMARKS 



0661 


51 


1507 


1550 


2004 


0662 


51 


1527 


1570 


2006 


0665 


54 


5000 


2100 


1650 


0664 


51 


2000 


2001 


2004 


0665 


50 


2002 


2100 


2006 


0666 


56 


2100 


2005 


2007 


0667 


54 


5000 


2100 


1607 


0670 


55 


0676 


2000 


0700 


0671 


55 


2000 


O676 


0675 


0672 


54 


0677 


2001 


0700 


0675 


56 


1550 


1657 


1550 


067U 


56 


1552 


1657 


1552 


0675 


54 


5000 


2100 


0605 


0676 


00 


0000 


0000 


0007 


0677 


00 


5200 


0000 


0650 


0700 


17 


2010 


5000 


0702 


0701 


54 


5000 


2100 


0706 


0702 


21 


1510 


2100 


0001 


0705 


21 


1551 


2100 


0001 


0704 


21 


1512 


2100 


0001 


0705 


21 


1555 


2100 


0001 


0706 


51 


1501 


1542 


1500 


0707 


51 


1521 


1562 


1520 


0710 


51 


1502 


15^5 


1501 


0711 


51 


1522 


1565 


1521 


0712 


51 


1510 


1551 


1502 


0715 


51 


1512 


1555 


1522 


0714 


51 


1550 


1571 


1505 


0715 


51 


1552 


1575 


1525 


0716 


51 


1506 


1547 


1505 


0717 


51 


1526 


1567 


1525 


0720 


51 


1507 


1550 


1506 


0721 


51 


1527 


1570 


1526 


0722 


51 


l 44 o 


l 4 oi 


1507 


0725 


51 


1442 


1405 


1527 


0724 


54 


5000 


2100 


0561 


0725 


55 


1512 


1455 


0751 


0726 


55 


0751 


1466 


0755 


0727 


52 


0750 


1677 


0757 


0750 


52 


0747 


l 64 o 


0740 


0751 


51 


0100 


0101 


2000 


0752 


25 


2001 


1551 


2005 


0755 


55 


2000 


1510 


2002 


0754 


51 


2002 


2005 


2002 


0755 


51 


0104 


0105 


2000 


0756 


54 


5000 


2100 


1540 


0757 


50 


2000 


2100 


0104 


0740 


50 


2001 


2100 


0105 


0741 


55 


0755 


1466 


0755 



Load 

f(2j) 

f(z 1+1 )/f(z i ) 

Load quotient 
for complex 
conjugate 
multiplication 

|f(z 1+ i)/f(z, )| ? 
is greater than 
100; halve 
X 1+ i &nd 
recompute z 1+1 
Floating point 
100 

If SEN 1 up, print iteration 

No, skip print 

Print 

z i+l 

real and 
imaginary 
Shift Zj..]_ 
to z ^„2 
Shift zi 
to ZJ..1 
Shift zi+i 
to Zj 

Shift X i+1 
to 

Shift f(zj„i) 
to ffZi.o) 

Shift f(z t ) 
to f(zi«l) 

Shift f(zi+i) 
to f(z^) 

Reiterate 

If iibag portion significant, 0751 
Initialize 0755 
Initialize 0757 
Initialize 0740 
Load A n 

} A n z i 

A n-1 + 

A n zi 

Put away 

terms of suppressed equation 
Modify to suppress 



IV - 1 - 8 



LAGRANGIAN INTERPOLATION METHOD FOP HIGH DEGREE POLYNOMIALS 



ADD 


INST 


M 1 


M 2 


M 3 


REMARKS 


0742 


55 


0737 


0264 


0737 


eqution to next 


07*0 


55 


0740 


0264 


0740 


lower degree 


0744 


54 


0656 


0735 


0732 


If extraction not complete, 0732 


0745 


36 


1372 


1331 


1372 


Reduce N* by 1 


0746 


34 


3000 


2100 


1013 


Transfer to conversion routine 


07b7 


00 


0000 


0000 


0103 


Constant for 0730 


0750 


00 


0000 


0000 


0102 


Constant for 0727 


0751 


35 


1310 


1331 


1303 


Build 2a term 


0752 


50 


1351 


2100 


1344 


in quadratic factor 


0755 


05 


3000 


3000 


1510 


Build 


0754 


31 


2002 


2003 


2006 


(a 2 + b 2 ) 


0755 


30 


2000 


2100 


2004 


term in 


0756 


36 


2100 


2001 


2005 


quadratic 


0757 


34 


3000 


2100 


1607 


factor 


0760 


31 


2000 


2001 


1304 




0761 


35 


0766 


1473 


2007 


Initialize 


0762 


55 


2007 


0264 


0770 


0770 


C765 


32 


1437 


1677 


1002 


Initialize 1002 


0764 


35 


1002 


1331 


2007 


Initialize 


0765 


32 


2007 


1535 


1003 


1003 


07 66 


31 


0100 


0101 


1300 


A n to TS 


0767 


31 


0102 


0103 


1301 


A n-1 4° b 


0770 


31 


0114 


0115 


1302 


A n-2 to TS c 


0771 


35 


1300 


1303 


2000 h 


I 


0772 


25 


1341 


1344 


2001 


x a to 200k , 2005 


0775 


31 


2000 


2001 


2004 J 


1 


0774 


35 


1300 


1304 


2000 


1 


0775 


25 


1341 


1345 


2001 


^ x ( a 2 + b 2 ) to 2006, 200? 


0776 


51 


2000 


2001 


2006 J 




0777 


31 


1301 


1342 


2000 


Load 


1000 


51 


1302 


1343 


2002 


^n-l 


1001 


34 


3000 


2100 


1410 


Add 


1002 


30 


1300 


2100 


0110 


Put away A n of 


1005 


30 


1341 


2100 


0111 


suppressed equation 


1004 


31 


2000 


2001 


1300 


Set up for 


1005 


31 


2002 


2003 


1301 


next term 


1006 


35 


0770 


1466 


0770 


Modify to pick up A n ^ 


1007 


35 


1002 


0264 


1002 


Modify put away 


1010 


35 


1003 


0264 


1003 


Modify put away 


1011 


34 


0656 


0770 


0770 


Has extraction been completed? 


1012 


36 


1372 


0264 


1372 


Yes, reduce N’ by 2 


1015 


30 


1351 


1310 


1301 


Form fractional part (l) 


1014 


56 


1510 


0261 


2000 


Build shift control word 


1015 


30 


1351 


2000 


2004 


Form integer part 


1016 


35 


2100 


2100 


1302 


Clear sign storage 


1017 


32 


1351 


II66 


1302 


Store sign 


1020 


34 


2004 


2100 


1102 


If integer ^ 0, 1102 


1021 


34 


3000 


2100 


1031 


No, 1031 



(l) The remainder of Channel 10 and Channel 11 is a binary floating 
point to decimal conversion which would be part of the overall program. 



IV - i - 9 



LAGRAHCIAN INTERPOLATION METHOD FOR HIGH DEGREE POLYNOMIALS 



ADD INST M 1 M 2 M 3 REMARKS 



1022 


12 


0000 


0000 


0000 


1025 


10 


0100 


0000 


0000 


1024 


00 


2400 


0000 


0000 


1025 


00 


1515 


2700 


0000 


1026 


10 


0000 


0100 


0000 


1027 


00 


3115 


1515 


1500 


1050 


30 


2005 


2100 


1300 


1051 


50 


1501 


2100 


2000 


1052 


30 


1167 


2100 


2007 


1053 


26 


2000 


0257 


2000 


1034 


27 


2007 


0256 


2007 


1035 


32 


2001 


0254 


2007 


1036 


37 


2007 


2100 


1035 


1037 


30 


2007 


2100 


1301 


lo4o 


52 


1016 


1022 


1301 


104 1 


34 


H63 


IO67 


1043 


1042 


54 


3000 


2100 


1170 


1043 


34 


1331 


2004 


1057 


1044 


27 


1101 


2100 


2007 


1045 


32 


1300 


0260 


2100 


1046 


34 


2000 


2100 


1053 


1047 


21 


1024 


1023 


0001 


1050 


30 


1300 


0256 


1300 


1051 


30 


2007 


0256 


2007 


1052 


54 


3000 


2100 


1045 


1055 


32 


1302 


1022 


1300 


1054 


21 


1300 


2007 


0001 


1055 


21 


1025 


1026 


0001 


1056 


34 


3000 


2100 


1066 


1057 


21 


1027 


1140 


0001 


1060 


32 


1302 


1166 


2100 


1061 


27 


2000 


0262 


2000 


1062 


34 


2000 


2100 


1065 


1063 


21 


1157 


1026 


0001 


io64 


34 


3000 


2100 


1066 


1065 


21 


ll60 


1026 


0001 


10 66 


21 


1301 


1101 


0001 


1067 


35 


1455 


1312 


1074 


1070 


21 


1161 


H65 


0002 


1071 


35 


1163 


2100 


IO67 


1072 


31 


1312 


1353 


1310 


1073 


34 


3000 


2100 


1013 


1074 


21 


1154 


1023 


0001 


1075 


34 


1372 


2100 


0310 


1076 


22 


0000 


0000 


0000 


1077 


35 


1164 


2100 


1067 


1100 


34 


3000 


2100 


1074 


1101 


01 


0000 


0000 


0001 



Suppress sign bits 
Print control 
Alphabetic print 
Alphabetic print 
Print control 
Alphabetic print 
Reset command 
Load fraction 
Load testword 

Multiply for most sig. digit 
Shift testword 

Extract digit into TW storage 
Repeat loop if o’flo 
Putavay converted fraction 
Suppress sign 
Print editing 
commands 
Print editing 
commands 

Extract most sig. dec. digit 
If 2000 ^ 0, 10^5 
If not print space 
Shift word 4 left 

and change print control 
Return to test next digit 
Enter sign 
Print integer 
Print decimal point 
Skip to 1066 
Print editing command 
Build sign print 

when only fractional 
part involved 
Print +. 

Skip to 1066 
Print 

Print fractional part 
If imag part insig, 1074 
No, print + j 
Set skip command in 1067 
Set up convert 
imaginary part 
Print carriage return 
If N * ^ 0, next root 
If N* * 0, halt 
Restore 1067 

Return to "test for end" 
Constant for 1044 



IV - 1 - 10 



LAGRANGIAN INTERPOLATION METHOD FOR HIGH DEGREE POLYNOMIALS 



ADD 


INST 


M 1 


M 2 


M 3 


REMARKS 


1102 


35 


2100 


2100 


2005 


High 


1103 


31 


1113 


2004 


2000 


integer 


llOl 


32 


2000 


11 14 


2100 


conversion 


1105 


30 


2000 


1115 


2000 


routine 


1106 


35 


1116 


2000 


1107 




1107 


35 


1126 


2100 


1110 




1110 


23 


2004 


11 40 


2000 




1111 


26 


2000 


1147 


2000 




1112 


54 


5000 


2100 


1142 




1115 


02 


0000 


0000 


0002 




1114 


00 


0000 


0000 


0174 




1115 


00 


0000 


0000 


0026 




1116 


35 


1115 


2100 


1110 




1117 


23 


2004 


1127 


2000 




1120 


23 


2004 


1130 


2000 




1121 


23 


2004 


1131 


2000 




1122 


23 


2004 


1132 


2000 




1125 


23 


2004 


1133 


2000 




1124 


23 


2004 


1134 


2000 




1125 


23 


2004 


1136 


2000 




1126 


23 


2004 


1137 


2000 




1127 


00 


0167 


1531 


1777 




1150 


00 


0013 


3274 


0777 




1151 


00 


0001 


1422 


2377 




1132 


00 


0000 


0750 


2177 




1153 


00 


0000 


0060 


6477 




1134 


00 


0000 


0004 


7037 




1135 


34 


3000 


2100 


1151 




1136 


00 


0000 


0000 


0307 




1137 


00 


0000 


0000 


0023 




11 4o 


10 


0000 


0000 


0100 




ll4i 


26 


2000 


1155 


2000 




1142 


30 


2005 


1150 


2005 




1143 


35 


1110 


1156 


1110 




1144 


35 


2001 


2005 


2005 




1145 


33 


1110 


1153 


ll4l 




1146 


34 


3000 


2100 


1030 




1147 


00 


0000 


0000 


0024 




1150 


00 


0000 


0000 


0004 




1151 


35 


2100 


2100 


2005 




1152 


34 


3000 


2100 


1146 




1153 


23 


2004 


1140 


2000 




1154 


00 


3400 


0000 


0000 




1155 


00 


0000 


OCOO 


0012 




1156 


00 


0000 


0001 


0000 





IV - 1 - 11 



LAGRANGIAN INTERPOLATION METHOD FOR HIGH DEGREE POLYNOMIALS 



ADD 


INST 


M 1 


M* 


M 3 


REMARKS 


1157 


00 


2352 


2764 


6400 


Print configuration for 1063 


1160 


00 


2252 


2764 


6400 


Print configuration for 1065 


1161 


00 


2515 


3621 


3264 


Print configuration for 1070 


1162 


00 


k26h 


6464 


3232 


Print configuration for 1070 


1163 


34 


3000 


2100 


1077 


Dummy reset for 1071 


ll64 


33 


1455 


1312 


1074 


Dummy reset for 1074 


1165 


10 


0000 


0000 


0000 


Print control - alphabetic 


1166 


02 


0000 


0000 


0000 


Sign extractor 


1167 


00 


0421 


0421 


0420 


Testword for 1072 


1170 


32 


1016 


1022 


1302 


Suppress integer sign when 


1171 


34 


2004 


2100 


1044 


printing imaginary portion 


1172 


21 


1174 


1175 


0001 


or print 0. 


1173 


34 


3000 


2100 


1066 


and return to 1066 


1174 


00 


5227 


0000 


0000 


Print configuration for 1172 


1175 


10 


0001 


0000 


0000 


Print control word 


1176 - 


• 1177 








Not required 


Channel 12 


is used 


for a memory check 


-sum routine in the root-solving 


routine and 


is available for other use when this program is incorporated 


in the 


overall synthesis program. 




Channel 13 


is used 


for temporary storage and for a few constants which 


are listed below. NOTE: This program 


is coded for minimum access only. 


i .e. , 


adjacent cells in the 


main memory are numbered 00, 4l, 02, 43, 


o4, 45 


, etc 


• 








1351, 


1557 


- floating point 


"1", must 


be addressed individually 


1533, 


1335 


- floating point 


"0", must 


be addressed individually 


1572 - 


N' - 


the degree of the reduced 


equation 


1374 


00 


3000 


3000 


0200 


Testword for 0505 


1576 


00 


0000 


0000 


0010 


Increment for 0303, 0304 


Channels l4 


, 15, and l6 contain complex floating point arithmetic 


routines, i 


.e., 










Complex Addition 




1410 




Complex Subtraction 




l4i6 




Simple Multiplication 


1500 




Simple Division 




1507 




Simple Subtraction 




1516 




Simple Addition 




1540 




Square Root 






1551 




Complex Multiplication 


1607 




Complex Division 




1650 



IV - 1 - 12 



LAGRANGIAN INTERPOLATION METHOD FOR HIGH DEGREE POLYNOMIALS 



ADD INST M 1 Vp REMARKS 

1400 - 1407 Used for temporary storage of operands 



l4io 


30 


3000 


1505 


1400 


l4ll 


32 


1400 


l64o 


1433 


1412 


04 


3000 


3000 


1400 


1413 


32 


1472 


1535 


1424 


l4i4 


32 


1471 


1677 


1430 


1415 


34 


3000 


2100 


1423 


l4i6 


30 


3000 


1603 


1536 


1417 


32 


1536 


1515 


1435 


1420 


04 


3000 


3000 


i4oo 


1421 


32 


1465 


I636 


1424 


1422 


32 


1467 


1636 


1430 


1423 


51 


2004 


2005 


2002 


1424 


34 


3000 


2100 


1516 


1425 


31 


2000 


2001 


1400 


1426 


31 


1402 


1443 


2000 


1427 


31 


l4o6 


1447 


2002 


1430 


34 


3000 


2100 


1516 


1431 


31 


2000 


2001 


2002 


1432 


31 


i4oo 


l44i 


2000 


1433 


34 


3000 


2100 


0624 


1434 - 


1436 


Temporary storage 


1437 


00 


0000 


0000 


0100 


1440 - 


1453 


Temporary storage 


1454 


02 


0000 


0000 


0044 


1455 - 


1463 


Temporary storage 


1464 


00 


0000 


7777 


7777 


1465 


00 


0000 


0000 


1516 


1466 


00 


0002 


0002 


0000 


1467 


00 


0000 


0000 


1516 


1470 


31 


0102 


0103 


2000 


1471 


00 


0000 


0000 


1540 


1472 


00 


0000 


0000 


1540 


1473 


00 


0004 


0004 


0000 


1474 


00 


0002 


0000 


0000 


1475 - 


1477 


Temporary storage 


1500 


30 


3000 


1505 


2006 


1501 


32 


2006 


1677 


1550 


1502 


35 


2000 


2004 


2000 


1503 


25 


2001 


2005 


2001 


1504 


34 


3000 


2100 


1547 


1505 


02 


0000 


0000 


0030 


1506 


00 


0000 


0000 


3777 


1507 


30 


5000 


1676 


2006 


1510 


32 


2006 


1677 


1550 


1511 


36 


2000 


2004 


2000 


1512 


23 


2001 


2005 


2001 


1513 


37 


2001 


3000 


1530 


1514 


3^ 


3000 


2100 


1547 


1515 


00 


0000 


0000 


7777 








IV - 


1 - 15 



Prepare 

exit 

Store operands 
Set 1424 

and 1430 for addition 
Transfer to 1423 
Prepare 
exit 

Store operands 
Set 1424 

and 1430 for subtraction 
Reposition real operands 
Transfer for ♦ reals 
Store sum of reals 
Load imaginary 
operands 

Transfer for + imaginaries 
Reposition imaginary sum 
Load real sum 
Exit 

Decrementer for 0306, 0307 



M 2 , M 3 extractor for 0333 
Constant for 1421 
Constant for 0323 
Constant for 1422 
Dummy reset for 0326 
Constant for l4l4 
Constant for 1413 
Constant for O332, 0336 
Incrementer (various) 

Prepare 

exit 

Add exponents 
Multiply magnitudes 
Transfer to exit 
Shift control for 1500 
M 3 extractor 
Prepare 
exit 

Subtract exponents 
Divide magnitudes 
If overflow, 1530 
Transfer to exot 
M 3 extractor 



LAGRANCIAN INTERPOLATION METHOD FOR HIGH DEGREE POLYNOMIALS 



ADD INST M 1 M 2 M 3 REMARKS 



1516 


30 


5000 


1505 


2006 


1517 


32 


2006 


1677 


1550 


1520 


36 


2100 


2003 


2003 


1521 


3k 


5000 


2100 


15 k2 


1522 


36 


2000 


2002 


2006 


1523 


30 


2001 


2006 


2007 


152 k 


30 


2002 


2100 


2000 


1525 


35 


2005 


200? 


2001 


1526 


57 


2001 


3000 


1550 


1527 


5k 


3000 


2100 


15k7 


1530 


27 


2001 


l6kk 


2001 


1531 


30 


2001 


l6kk 


2001 


1532 


27 


2001 


l6k6 


2001 


1533 


35 


2000 


l6k6 


2000 


153 k 


3k 


3000 


2100 


15k7 


1535 


00 


0000 


0000 


7777 


1536 


00 


0000 


0000 


062 k 


1537 


02 


0000 


0000 


0001 


15^0 


30 


3000 


1505 


2006 


15 kl 


32 


2006 


1515 


1550 


15 k2 


33 


2002 


2000 


1522 


15^3 


36 


2002 


2000 


2006 


15 kk 


30 


2003 


2006 


2007 


15^5 


35 


2001 


2007 


2001 


15 k 6 


37 


2001 


3000 


1530 


15^7 


31 


2000 


2001 


2000 


1550 


3k 


3000 


2100 


• • « • 


1551 


50 


3000 


1676 


2006 


1552 


52 


2006 


1535 


1550 


1553 


35 


2000 


0261 


2000 


155k 


32 


2000 


1331 


210k 


1555 


30 


2000 


l6kk 


2006 


1556 


3k 


200k 


2100 


1561 


1557 


30 


2001 


l6kk 


2007 


1560 


3k 


3000 


2100 


1562 


1561 


30 


2001 


2100 


2007 


1562 


30 


1575 


2100 


2002 


1565 


23 


2007 


2002 


2001 


156 k 


36 


2001 


2002 


2001 


1565 


30 


2001 


l6kk 


2001 


1566 


3k 


2001 


2100 


1570 


1567 


3k 


3000 


2100 


1572 


1570 


35 


2002 


2001 


2002 


1571 


3k 


3000 


2100 


1563 


1572 


25 


2002 


1577 


2001 


1573 


36 


2006 


1576 


2000 


157k 


3k 


3000 


2100 


15 k7 



Prepare 

exit 

Change sign 2 nd operand 
Transfer into add 
Form shift control 
Shift 1®^ operand 
Shift 2 exponent to Ans. cell 
Add operands 
If overflow, 1530 
Transfer to exit 
Normal 
correct 
overflow 
procedure 
Transfer to exit 
M 3 extractor 
Prepare exit storage 
Shift control for 1555 
Prepare 
exit 

If 2002 ^ 2000, 1522 

No, form shift control 

Shift 2 nd operand 

Add operands 

If overflow, 1530 

Normalize 

Exit 

Prepare exit 

from square root routine 
Multiply by 2 44 

Extract last bit of exp into 200k 
Take square root of exp 
If exp odd, magnitude has been 
divided by 2, by 1555 
If exp even, divide magnitude 
by 2 

Skip to iteration 

Load magnitude for iteration 

Load a s 

(x/2)/a. 

ix/2ei/) - aj^ 

l/2(x/2a 1 ) - ai = a i+ i - ^ 

Has process converged? 

Yes, 1572; no, 1570 

a i+l a i 

and reiterate 

Multilpy by 42/2 
Divide by 2 21 
Transfer to exit 



IV - 1 - lk 



LAGRANGIAN INTERPOLATION METHOD FOR HIGH DEGREE POLYNOMIALS 
ADD INST M 1 M 2 M 3 REMARKS 



1575 


00 


7777 


7777 


7777 


Constant for 1562 


1576 


00 


0000 


0000 


0021 


Constant for 1573 


1577 


00 


5520 


2363 


1500 


Constant for 1572 



1600 - 1602 Used for temporary storage of operands 

1603 02 0000 0000 0030 Shift control for 1607 

1604 - 1606 Used for temporary storage of operands 

Prepare 



1607 


30 


3000 


1603 


1601 


1610 


52 


1601 


1636 


1635 


l6ll 


04 


3000 


3000 


1600 


1612 


5^ 


5000 


2100 


1500 


1613 


30 


2000 


2100 


1601 


l6lU 


30 


2001 


2100 


1605 


1615 


35 


1602 


1606 


2000 


1616 


25 


1643 


1647 


2001 


1617 


51 


2000 


2001 


2002 


1620 


51 


1601 


1605 


2000 


1621 


5^ 


3000 


2100 


1516 


1622 


30 


2000 


2100 


1601 


1623 


50 


2001 


2100 


1605 


1621* 


35 


1602 


1604 


2000 


1625 


25 


1643 


1645 


2001 


1626 


31 


2000 


2001 


2002 


1627 


55 


1600 


1606 


2000 


1630 


25 


l64i 


1647 


2001 


1631 


51 


2000 


2001 


2000 


1632 


3^ 


3000 


2100 


1540 


1633 


31 


2000 


2001 


2002 


1634 


31 


1601 


1605 


2000 


1635 


5^ 


3000 


2100 


0627 


1636 


00 


0000 


0000 


7777 


1637 


00 


0000 


0000 


0001 


l64o 


00 


0000 


0000 


7777 


l64i 


Used 


for temporary 


storage 


1642 


00 


0000 


0000 


0001 


1643 


Used 


for temporary storage 


1644 


02 


0000 


0000 


0001 


1645 


Used 


for temporary storage 


1646 


00 


0000 


0000 


0001 


1647 


Used 


for temporary storage 


1650 


30 


5000 


1676 


1601 


1651 


32 


1601 


1677 


1635 


1652 


o4 


3000 


3000 


1600 


1653 


30 


2004 


1637 


2000 


1654 


25 


2005 


1645 


2001 


1655 


31 


2000 


2001 


2000 



exit 

Store operands 
Multiply reals 
and store results 
in 1601, 1605 
Multiply lmaglnarles 
to obtain negative real 
and shift for subtraction 
from previous real product 
Transfer to subtract 
Put away real product 
in 1601, 1605 
Multiply for 

first cross-product tern 
and hold in 2002, 2003 
Multiply second 

cross-product term 
and prepare for addition 
Add for imaginary cross-products 
Put results in2002, 2003 
Insert real product 
Exit 

M 3 extractor for l6l0 
Shifter for 1656 
M 3 extractor 

Shift control word 

Shift control word 

Shift Control word 

Prepare (A +JB) 

exit (C +JD) 

Store operands 
? orra 
C 2 

and normalise 



TV - 1 - 15 



LA GRAN GIAN INTERPOLATION METHOD 



ADD 


INST 


M 1 


M 2 


M 3 


1656 


50 


2006 


1637 


2002 


1657 


25 


2007 


1647 


2003 


1660 


31 


2002 


2003 


2002 


1661 


34 


3000 


2100 


1540 


1662 


31 


2000 


2001 


2004 


1663 


51 


1604 


1645 


2000 


1664 


34 


3000 


2100 


1507 


1665 


30 


2000 


2100 


1604 


1 666 


30 


2001 


2100 


1645 


1667 


31 


1606 


1647 


2000 


1670 


34 


3000 


2100 


1507 


1671 


30 


2000 


2100 


1606 


1672 


36 


2100 


2001 


1647 


1673 


31 


1604 


1645 


2004 


1674 


31 


1600 


l64i 


2000 


1675 


34 


3000 


2100 


1612 


1676 


02 


0000 


0000 


0030 


1677 


00 


0000 


0000 


7777 



FINIS 



iv - l - 16 



FOR HIGH DEGREE POLYNOMIALS 
REMARKS 



Form 

D 2 

and normalize 
(C 2 + D 2 ) 

Shift to 2004, 2005 
Load C 

C/ (C 2 + D 2 ) 

Store in 
l6o4, 1645 
Load D 
D/ (C 2 + D 2 ) 

Store in 1606, 

1647 vrtth negative sign 
Load conjugate 
quotient and 
transfer to multiply 
Shift control vord 
M 3 extractor 



FLOW CHART OF LAG RANG IAN INTERPOLATION METHOD 
FOR SOLUTION OF ALGEBRAIC EQUATIONS OF HIGH DEGREE 







A 



\t + 1 — ^ 






> — 









IV - 2 - 1 















APPENDIX V 



USE OF THE EXISTING PROGRAMS ON THE CRC-102A 

All existing programs have been coded for minimum -access word channel 
and will not work in the sequential channels, 

1. Dirichlet Power Series Expansion 

Since the synthetic division and conversion program has not been 
written as yet, if f^(t) and f^(t) are given, ^h^ ^ must be computed by 
hand. If h(t) is given, it must be converted to ^h^ by multiplying by 
At. As the program is now located, storage is available for 100 (octal) 
ordinates in floating point form if the decimal to floating point con- 
version in Channel 2 is to be used. If the conversion is accomplished 
prior to entry of data, 150 (octal) ordinates may be used. 

T 

S = 2 — = 2 NT where S * total cells needed 

At 

T » number of time units problem is 
scaled over (note command 0423) 

At a sampling interval 

N a samples per time unit 

The program may be relocated to accommodate more sampled ordinates. 
The binary floating point numbers representing the ordinates of h(t) 
are entered in cells 0000 - 0117, exponents in even cells, magnitudes in 
odd; this allows 40 ordinates to be used. 

The test switches determine the amount of interpolation performed. 

The output will be H*(s) in decimal form for as many terras requested. 

The number of terms is set by entering this datum as an octal number in 
J, cell 0540. 

2. Matrix Inversion 

The matrix inversion program takes the ascending power H*(s) terms in 
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